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Emphasis on Heat-Sustaining Nose Shapes at 
Hypersonic Speeds' 


A. J. HANAWALT,* A. H. BLESSING,** anp C. M. SCHMIDT*** 
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SUMMARY 


The leading edges and noses of hypersonic vehicles are sub- 
jected to severe aerodynamic heating and must be cooled in some 
manner—e.g., internal convection, transpiration, or radiation. 
It is this latter mode of handling the problem that is discussed 
in this paper. Neglecting conduction in the leading-edge 
region, the maximum temperature for long-range hypersonic 
gliders is of the same order as the melting point of refractory 
materials, with a corresponding large temperature gradient 
away from the leading edge. Inclusion of conduction in the 
aft direction reduces the maximum temperature and distributes 
the heat to a location that will radiate it out from the surface. 
For either steady-state or transient conditions, the temperature 
at the leading edge is reduced by conduction, while the tempera- 
ture aft of the leading-edge shoulder is increased, thus setting 
up a heat transmission balance between the convective influx 
of heat, the redistribution of heat by conduction, and the radia- 
tion of heat from the surface. The feasibility of such a mecha- 
nism can be enhanced by suitably choosing leading-edge 
shapes and materials. The philosophy behind the choice of 
leading-edge shapes is discussed and the effects of varying 
parameters, such as shape, diameter, emissivity, conductivity, 
thickness, etc., are shown. 


SYMBOLS 


lift coefficient 
diameter, ft. 

= gravitational constant, ft./sec.? 
enthalpy, B.t.u./Ib. 
heat-transfer coefficient, B.t.u./ft.2hr.°F. 
thermal conductivity, B.t.u.-in./ft.? hr.°F. 
Mach Number 
coordinate normal to the outer surface 
Prandtl Number 
pressure, Ib. /ft.? 


+ This paper is based on a paper presented at the 1958 Heat 
Transfer & Fluid Mechanics Institute, June 19, 20, and 21, 
1958. Received July 8, 1958. Revised and received October 
23, 1958. 

* Head, Heat Transfer Group, Space Flight Division. 

** Research Engineer, Heat Transfer Group. 

*** Systems Engineer. 


= convective heat flux, B.t.u./ft.*hr. 
local body radius, ft. 
radius, ft. or gas constant, 53.3 ft.-lb./Ib.°F. 
reference area, 
arc length, ft. 
temperature, °F. or °R. 
recovery temperature, °R. 
velocity, ft./sec. 
satellite velocity, ft./sec. 
weight of vehicle, Ib. 
Cartesian coordinates, x aligned with body axis, ft. 
distance from orbit to center of earth, ft. 
semi-wedge angle, deg. 
emissivity 
sweep angle, deg. 
viscosity, slugs/ft. sec. 
density, slugs/ft.* 
Stefan-Boltzmann constant, 1.73 B.t.u./ 

ft.2hr.(°R.)4 
T = skin thickness, ft. 


one 


oF 


Subscripts 
local conditions outside the boundary layer 
stagnation condition 
stagnation conditions behind the normal shock 
evaluated at the wall temperature 
free-stream conditions 


(1) INTRODUCTION 


- at velocities many times the speed of sound 
has become a matter of considerable concern. One 
of the attendant problems of such high-speed flight 
is aerodynamic heating. Its severity has been ade- 
quately demonstrated in the literature, and, in fact, 
for most advanced aircraft it is one of the dominant 
design considerations. It is the purpose of this paper 
to treat one major aspect of this problem, suggest a 
solution, demonstrate trends, and indicate various 
important parameters in this type of solution. 

The aerodynamic heating problem can be met in 
several ways. Two different methods of solution 
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have been attempted for the long-range ballistic missile 
with its formidable re-entry heating problem—.e., 
heat sink and ablation cooling. Both of these methods, 
particularly ablation, are well adapted to short periods 
of very high heat flux, but would not be suited for 
long flight times. Transpiration and film cooling are 
theoretically possible for flights of both short and long 
duration but have two principal disadvantages—(1) a 
materials problem and (2) a strong tendency to cause 
transition from laminar to turbulent flow. A fourth 
method of cooling, radiation, can be accomplished by 
allowing the surface temperature to reach equilibrium— 
i.e., the temperature at which the emitted radiation 
just balances the convective influx of heat. The glide 
or cruise vehicle falls into the category of vehicles 
that can be cooled in this manner. 

The intention of this paper is to discuss the heating 
problem of a long-range hypersonic glider, with atten- 
tion focused on the maximum heat input region—.e., 
the entry edges. It has been adequately shown that, 
for maximum glider range, the lift to drag ratio should 
be a maximum, which in turn means that the leading 
edges and nose should be small to reduce drag. This 
means high heating rates and radiation equilibrium 
temperatures from 4,000 to 5,000°F., depending on the 
range of the glider, its wing loading, and its lift coeffi- 
cient. In general, these temperatures are beyond the 
limits of presently available structural materials and 
require some means of reduction. Insulation or 
ablation will not work because of the long time of 
exposure to such high temperatures—half an hour or 
more. The use of transpiration or film cooling is 
objectionable because of the limitations of materials 
and the importance of maintaining laminar flow 
further back on the glider. To simply add enough 
mass so as to absorb all the incoming heat would 
require large leading-edge diameters—i.e., high drag— 
and excessive weight. The use of an internal coolant, 
either forced convection or pool type boiling, shows 
considerable promise and would be feasible; however, 
this approach is neither as simple nor as light as 
desired. 

A more detailed examination of radiation cooling, 
including the mechanism of conduction, shows that 
the heat-sustaining leading edge is more nearly prac- 
tical than at first suspected. In fact, for reasonable 
values of conductivity, the maximum equilibrium 
temperatures can be reduced to an acceptable level for 
presently available heat-sustaining materials. It is 
the purpose of the present paper to investigate the 
degree of reduction that can be obtained plus the rela- 
tive importance of the various parameters. This has 
been done for the maximum heating condition of a 
glider as determined in Section (2). The prediction 
of aerodynamic heating rates is treated in Section (3), 
where the methods of analysis are discussed and 
referenced. Section (4) is concerned with the problems 
of obtaining numerical results, either by hand or high 
speed computers, while Section (5) presents the results 
of the study and investigates the effects of varying 
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such parameters as conductivity, flight altitude, wing 
sweep, and diameter. 


(2) Maxtmum HEATING CONDITION FOR A HYPERSONIC 
GLIDER 


The following is an approximate derivation of the 
location in the equilibrium glide trajectory where 
maximum heating conditions are encountered. For 
a glider which flies an equilibrium path, where at all 
times lift equals apparent weight, maximum heating 
does not occur at maximum velocity. The higher 
the velocity the greater the centrifugal force and there- 
fore the smaller the apparent weight, or lift. Since 
aerodynamic heating is related to lift, it can be shown 
that the maximum heating rate occurs at some velocity 
less than orbital. This is demonstrated by the follow- 
ing equations, where Eq. (3) presents the velocity at 
which maximum heat flux occurs: 


q = h(T, — Ts) « VPV V? (1) 
P « Lift « (W/S)[1 — (V?/V,?)] (2) 
V = V5/7 V, & 22,000 ft./sec. (3) 


The corresponding altitude can be obtained from the 
following equations and is seen to be a function only 
of wing loading and lift coefficient: 


Centrifugal Force + Lift = Weight (4) 
W(V2/g20) + = W (5) 
p = — (1/gz0)] (6) 


For a C, = 0.5 and W/S = 50 lb./ft.?, the equilibrium 
altitude is 241,000 ft. This combination of velocity 
and altitude serves as a basis for computation in the 
following Sections. Any other point in the glide 
trajectory would present a less severe heating condition. 


(3) CONVECTIVE HEATING RATES 


Before a heat-sustaining leading edge can be de- 
signed, it is necessary to describe accurately the 
thermal environment and the air loads to which the 
leading-edge region is subjected. Of these, aero- 
dynamic heating is the most critical and one which 
has received considerable attention in the literature. 
The earliest attempts at analysis considered the flat 
plate in incompressible flow. Later this was extended 
to the stagnation points of cylinders and spheres by 
Squire’ and Sibulkin,‘ respectively. More recently, 
the stagnation point has been analyzed including the 
effects of compressibility? and even real gas effects 
such as dissociation.’ ® Needless to say the latter 
three methods are somewhat similar and roughly 
predict the same results. Because of the fact that 
reference 5 applies to all flight speeds and wall tem- 
peratures with no additional modifications, it is the 
method used for this analysis. The other methods 


would give essentially the same results and permit 


the same conclusions. 
The distribution of the convective heat flux imme- 
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diately aft of the stagnation point is less adequately 
understood. Two of the methods that come closest 
to describing the physical situation are those of Lees! 
and Stine and Wanlass.’ The former neglects the 
effect of the pressure gradient tern in the boundary- 
layer equations, while the latter assumes local siimi- 
larity. Because of its greater simplicity and also 
that it generally gives similar results, Lees’ method 
was used to deterinine the ratio of ¢/Qsiay aft of the stag- 
nation point with reference 5 used to compute Qgrag. 
Using the relation pu = constant across the boundary 
layer, Lees’ method reduces to 


= 0.51Pr7?? (Ho — He) 


where m = O for two-dimensional flow and m = | for 
axisyininetric flow. The most deficient part of this 
method is its presumption of a knowledge of the flow 
conditions, since neither theory nor experiment exists 
which accurately provides the necessary values. An 
approximate imethod of accounting for inviscid blunt 
nose effects, which predicts much higher pressures aft 
of the shoulder than does Newtonian flow, was used to 
predict the pressure distribution. This approximate 
method was obtained using the theories of references 
8 and 9. 

Application of Eq. (7) to a hemisphere-cylinder 
produced the heat flux distributions shown in Fig. 
l(a). It is seen that ¢/Qstay is a strong function of 
Mach Number with the variations in this ratio being 
greatest for the higher Mach Numbers. Eq. (7) was 
prograined for an IBM 704 from which Fig. 1(a) was 
obtained. Similar curves were obtained for a hemi- 
cylinder-slab and were used for the following computa- 
tion of temperatures. Fig. 1(a) shows that, regardless 
of size, the heating rates over a hemisphere vary by 
an order of magnitude from stagnation point to 
Shoulder. Fig. 1(a) also shows that ¢/qs:ay is a function 
only of the ratio of the arc length to diameter or radius. 
Since gs1a9, and therefore g, is inversely proportional to 

D, decreasing the diameter increases the value of 
dq/ds by the factor 1/D*’?. Thus, for small diameters 
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dq/ds is very large and, if conduction is neglected, 
large variations in the radiation equilibrium tem- 
perature will result—as much as 1,500°F. in less than 
an inch. This, of course, suggests the fallacy of 
neglecting conduction from which large benefits can 
be expected to be obtained. 

One way in which the importance of conduction can 
be accentuated is shown in Fig. 1(b), which presents 
the heat flux distribution for a hemicylinder-slab and 
a “blunt” leading edge. In order to show the proper 
relationships, gsi, for the “‘blunt’’ leading edge has 
arbitrarily been taken to be the heat flux at the stag- 
nation point of the corresponding diaimeter cylinder. 
It is seen from Fig. 1(b) that, by proper choice of shape, 
the maximum heating rate is not altered much; how- 
ever, it is located near a region of low heat input, thus 
shortening the conduction path. The benefits to be 
gained from the use of such shapes are exemplified 
below for both the two-diinensional and three-dimen- 
sional bodies. 


(4) Heat CONDUCTION IN THE LEADING EDGE 


With the foregoing expressions for the heat flow to 
a leading-edge region, a semiparametric investigation 
of the temperature distribution in a practical leading 
edge was undertaken. Since it was deemed advisable 
to determine maximum temperatures in this region, 
equilibrium (steady-state) conditions were examined. 
This is particularly appropriate because of the stated 
application for a hypersonic glide vehicle for which 
flight conditions change slowly. In order to accomplish 
this investigation, the Fourier heat conduction equation 
with the following boundary conditions must be 
solved: 


+ (0?7/dy?) = 0; 
—k(OT/On) = g — eoT ‘at the external surface; (8) 
OT /Ox = 0 at the rear of the leading-edge region, i 


If the leading edge is at a zero angle of attack, the heat 
flux into the upper and lower surface will be equal; 


hence, the temperatures will be symmetric with 
1.0 ] ] ] 
\ | 
CYLINDER 
6 A 
wa | 
BLUNT LE. - 
4 
| 
O 2 4 6 8 1.0 
s/D 
Fic. 1(b). Heat flux distribution on a leading edge. 
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respect to the centerline. For this condition only half 
of the leading edge need be examined; this is done by 
assuming that the centerline is insulated. Eq. (8) 
js not amenable to analytic solution because of the 
complex form of g (a function of diameter, shape, sweep, 
wedge angle, and distance aft) and the fourth degree 
radiation term. As a result, numerical techniques 
must be considered. The particular one employed 
here is the simple finite difference method, which can 
be readily programed for high-speed computing ma- 
chines so that a large parametric investigation can be 
done very quickly. It consists of subdividing the 
leading edge into N elements, not necessarily of uni- 
form size. Eq. (8) is then replaced by a heat balance 
equation for each element, involving both convection 
and radiation terms. Because of the fourth degree 
temperature term these equations cannot be solved 
directly; however, the solution can be obtained by an 
iterative process if the fourth degree terims are approxi- 
mated by the first two terms of a Taylor series ex- 
pansion—.e., 


where 7';* is a guessed temperature for the ith element. 
It was found during the course of this investigation 
that it was not difficult to guess 7;* within 200°F. of 
the correct 7;, which meant an error of less than 
3 per cent in 7;*. Using this approximation the 
resulting V simultaneous linear equations in as many 
unknowns can be solved simply and quickly by a 
matrix inversion routine on a digital computer, such 
as the IBM 704. The only time-consuming portion 
of this analysis is the determination of the coefficients 
that are supplied to the machine. However, it would 
be entirely feasible to write a program so that the 
machine can compute these coefficients itself. Until 
some experience is obtained, it is suggested that the 
solution resulting from the first guess be used as a 
second guess and the procedure repeated. It was 
found that this process is rapidly convergent and that 
in most cases a second iteration was not required. 
This depends, of course, on being able to make a good 
first guess and on the accuracy required. 

If a high-speed computer is not available, results 
can be obtained with a slide rule, for one-dimensional 
conduction, by a different process. The leading edge 
is again broken up into elements, heat balance equa- 
tions are written for each element, the temperature of 
the first element is guessed, and the temperature of 
the second element is computed from a linear equation. 
This is continued for successive elements until the 
aft boundary condition is obtained. The process is 
repeated until the correct boundary condition is 
satisfied. 

For most of the cases examined, a one-dimensional 
analysis was utilized because the skin thicknesses 
were sufficiently small that rather insignificant gradients 
existed across the wall. This approach gives results 
that are very representative of the average tempera- 
ture through the thickness. For a specific case, a 
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two-dimensional analysis was made in order to de- 
termine heat flux lines and isotherms. Except for a 
very small region near the stagnation point the results 
agreed very well with the one-dimensional analysis. 


(5) RESULTS 


As mentioned earlier, a semiparametric investigation 
was performed to reveal trends and important variables. 
The base point conditions for this analysis were: 


Velocity, 22,000 ft. sec. 

Altitude, 241,000 ft. 

Lift coefficient, 0.5 

Shape, Circular cylinder followed by a solid slab wing 
Radius, 0.5 in. 

Conductivity, 250 B.t.u.-in. ft.*hr.°F. 

Sweep angle, 0° 

Emissivity, 0.9 

Length of leading-edge region, 6 in. 


Because of the large number of parameters involved, 
the investigation was limited to only a few values of 
each and then only as branches from the base point. 
For the conditions considered, the maximum tempera- 
ture was rather insensitive to the length of structure 
perpendicular to the leading edge and to the boundary 
condition at the aft end—e.g., either insulated or the 
no-conduction equilibrium temperature. For the most 
part, the present study assumed that the structure was 
insulated at the aft end of the leading-edge region. 

The results of the analysis for the base point are 
shown in Fig. 2 as temperature versus distance aft, 
along with the results for several different values of 
conductivity. The benefits to be gained by including 
conduction are obvious, not only for reducing the 
maximum temperature but also the temperature 
gradients. The k = O curve implies no-conduction 
and is seen to have not only the highest temperature 
but also the largest temperature gradient. For k = 
250 B.t.u.-in. /ft.*hr.°F., representative of graphite 
or silicon carbide at these temperatures, 7), is reduced 
920°F.; in addition, the temperature gradients are 
much smaller. The k = © curve means an isothermal 
body and represents the largest reductions in Tnoz 
attainable by conduction. Unlike that for the k = 0 
or 250 cases, the actual temperature for k = © isa 
strong function of just how far aft conduction takes 
place. Fig. 3 is a cross plot of the maximum tempera- 
tures shown in Fig. 2 and reveals a rather sharp knee 
in the curve at about k = 500 B.t.u.-in. ft.*hr.°F. 

The effect of lift coefficient on the maximum tem- 
perature is demonstrated in Fig. 4. As shown in 
Eq. (6), flight altitude is related to lift coefficient. 
For example, increasing C, from 0.5 to 1.0 increases 
the flight altitude from 241,000 to 255,000 ft. This 
means a decrease in density and, therefore, in heat- 
transfer coefficient, which not only decreases 7 for 
no-conduction but also when conduction is included. 
The important result shown in Fig. 4 is that the no- 
conduction and k = 250 curves parallel each other 
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very well, indicating that the significance of conduction 
is rather invariant with changes in flight altitude— 
i.e., convective heat-transfer rate. The effect of 
wedge angle or angle of attack variation is shown in 
Fig. 5 and can be seen to be virtually nonexistent. 
This is because the heating rate on the wing is an order 
of magnitude less than it is at the stagnation point, 
and, hence, even if the flat-plate heating rate is doubled, 
the maximum temperature is affected only slightly. 
It should be noted that Fig. 5 was prepared for the 
case of a given skin thickness; if the wing were kept 
solid, the variation would be even less. 

Sweeping the wing leading edge is a well-known 
method of reducing the heating rates to and the 
temperatures of leading edges. It has been shown 
that the heat flux varies essentially as the cosine of the 
angle of sweep from 0 to 60°. The effect that this 
has on the maximum equilibrium temperature includ- 
ing conduction and no-conduction is demonstrated 
in Fig. 6. It can be noted that sweeping the wing up 
to 30° is not very helpful in reducing the temperature, 
but beyond this point it becomes increasingly ad- 
vantageous. Again the curves are very nearly parallel 
so that the conclusions that are presented for A = 0 
are equally applicable for other sweep angles as long 
as the temperature level is reduced accordingly. 

The variation of maximum temperature with skin 
thickness is shown in Fig. 7 and can be seen to be 
rather linear. Zero thickness implies no-conduction 
while 0.5 in. thick means a solid wing, a thickness 
which may not be exceeded as long as this particular 
diameter leading edge is being examined. 

The effect of leading-edge diameter is demonstrated 
in Fig. 8. Because the heat flux rate is proportional 
to the inverse of the square root of the diameter, the 
no-conduction equilibrium temperature increases with 
decreasing diameter and would tend to approach the 
stagnation temperature as the diameter goes to zero. 
However, the maximum temperature with conduction 
behaves very differently, as is seen in this Figure. 
For large diameters, the conduction path is too long 
to offer any significant reduction from the no-conduc- 
tion case. As the diameter becomes smaller, con- 
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duction plays an ever increasingly important role in 
reducing the temperature because of the shorter 
conduction paths; in fact, a slight inversion in the 
Tmax Curve takes place at a diameter of about two in. 
For infinite conductivity, this inversion would occur 
at a much larger diameter. 

Up to this point, only circular cylinder leading edges 
have been considered. However, blunter leading edges 
have shown a good deal of promise in flight tests and 
from limited analytical work. Not only do the 
blunter shapes show benefits from the standpoint of 
less total heat flux, but also the conduction path is 
enhanced a good deal. The maximum heating rate 
on a hemicylinder occurs at the stagnation point; 
however, for blunter shapes this peak value moves 
toward the shoulder. The drop in heating rate from 
the maximum heating point back along the afterbody 
is, therefore, much more pronounced than on the 
hemicylinder-slab, and, thus, the potential for con- 
ducting heat away is increased. An example of the 
benefits to be gained is shown in Fig. 9, where a 1/10 
power leading edge is compared with a hemicylinder. 
The additional reduction in 7); is large enough to 
show that this is the direction to go with leading-edge 
shapes. 

Temperatures on the fuselage nose will tend to be 
higher because of two factors—higher heating rates 
to the corresponding axisymmetric body and _ the 
unavailability of the alleviating sweep effect. On 
the other hand, relatively large noses do not contribute 
excessively to the drag of a glider, which means that 
large diameters could be employed to minimize the 
temperatures. Since this is the case, appreciable 
benefits can be obtained with conduction and a proper 
choice of nose shape. Fig. 10 shows predicted tem- 
peratures for three nose shapes with a 1.0 ft. diameter. 
Although the only effect of conduction for the hemi- 
sphere-cone is to produce a slight smoothing of the 
temperature gradient at the shoulder, it does reduce 
the maximum temperature for the other two shapes 
by 275°F. and 360°F., respectively. It would be 
possible to obtain even larger reductions from more 
optimum nose shapes. 
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Fic.9. Effect of leading-edge shape on temperature distribution. 
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THERMAL ANALYSIS OF 


The surface emissivity is of course a very important 
parameter. For values less than 0.9, the maximum 
temperature is, of course, greater. However, the 
difference between 7,,,., for conduction and no-conduc- 
tion increases with decrease in emissivity, indicating 
the greater importance of conduction for the lower 
surface emissivities. 

For materials like graphite, the thermal conductivity 
is a strong function of temperature. This means that 
Eq. (8) should be replaced by the more general form 
of the Fourier heat equation: 


(0°?7 Ox?) + (0°?T/Oy?) + 
(1/k)(Ok/OT) [(OT/Ox)? + (OT /dy)?] = (10) 


A sample calculation was performed including the 
term Ok/OT for graphite, and it was found that Tino: 
decreased very slightly, less than 50°F. 

The present analysis has concerned itself only with 
heating rates and temperatures. This was done for 
the purpose of showing the importance of conduction 
as a means of reducing the temperature level. No 
attention was given to thermal stresses, weight, 
stability, and control as affected by wing sweep, or 
changes in flight attitude necessary to effect changes 
in lift coefficient. In an actual design, it would be 
necessary to consider all of these in order to better 
design the system. 


(6) CoNcLUSIONS 


The leading edge of a hypersonic glider was examined 
for maximum heating rates and temperatures. The 
importance of conduction was demonstrated by 
numerical examples, which showed that, if graphite 
could be used, the maximum temperature would be 
reduced by 900°F. from the no-conduction value of 
4,300°F. Trends and important variables were re- 
vealed by a semiparametric analysis, which investigated 
such things as conductivity, lift coefficient, skin thick- 
ness, wedge angle, diameter, sweep angle, and shape. 
Of these, only wedge angle and diameter were found 
to be unimportant. Conductivity was shown to be 
the most important single variable with large benefits 
being gained up to a k of 500 B.t.u.-in./ft*hr.°F. 
Large lift coefficients and sweep angles were found to 
be beneficial in reducing the actual magnitude of the 
temperature, but, because the no-conduction tempera- 
tures are also affected, the importance of conduction 
remained relatively the same. Shape was also revealed 
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to be an important variable because of its effect on the 
location of the maximum heating point and the result- 
ing reduction in conduction path. In particular, a 
1/10 power leading edge showed a 250°F. gain over a 
hemicylinder, with similar gains being obtained for 
a large diameter nose. 
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The Characteristics of Roughness From 
Insects as Observed for Two-Dimensional, 
Incompressible Flow Past Airfoils' 


WALTER S. COLEMAN* 
Blackburn and General Aircraft, Lid. 


SUMMARY 


Advances in the practical development of boundary-layer 
control for the maintenance of extensive laminar flow have drawn 
attention to the problem of surface roughness, due not only to 
artificial irregularities such as rivet heads, lap joints, window 
panels, etc., but also to the kind generated in flight from impact 
with insects. This natural form of roughening, the effects of which 
have been noted, though not investigated previously, is the sub- 
ject of the present paper. 

The phenomenon may be divided into two parts—namely, 
(1) the manner of distribution and magnitude of the roughness, 
and (2) its effect upon the stability of the laminar boundary 
layer. Wind-tunnel experiments with the fruit fly, Drosophila, 
and the common housefly for the investigation of both (1) and 
(2) in the case of two-dimensional, incompressible flow past 
airfoils are fully described. The former problem has also been 
treated mathematically in a separate paper, not yet published, 
agreement between theory and experiment being satisfactory in 
all essentials. 

The characteristics of the roughness profile consist principally 
of a pronounced peak near the leading edge, followed by an ex- 
tensive area of surface over which there is a much reduced and 
gradually diminishing value of the excrescence height. Further, 
it is shown that, if the severe leading-edge roughness, or its effect 
upon the boundary layer, can be eliminated, then the down- 
stream roughness causes no disturbance to the passage of a lami- 
nar layer—i.e., the surface, though roughened, is aerodynamically 
smooth. Moreover, it appears that the conditions defining the 
upstream boundary to this region of insignificant roughness are 
fundamentally the same as those which determine the critical 
state for transition at an artificial disturbance of a three-dimen- 
sional character. 


SYMBOLS 


x = coordinate along the chord from the leading edge 

s = coordinate along the surface from the leading 
edge 

ee = total extent of the roughness measured along the 
chord and surface, respectively 

is Se = position of transition measured along the chord 


and surface, respectively. Also equivalent to 
the significant extent of natural roughness 
[see Section (6)] 


Ax = amount of natural roughness removed from the 
surface [see Section (6)] 

n = coordinate normal to the surface, measured 
positively into the fluid from origin on the 
surface 
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c = chord of airfoil 

€ = excrescence height at any point on the surface, 
taken as the value defined by the envelope 
[see Section (4)] in the case of natural rough- 
ness 

= critical excrescence height for premature transi- 


Ccrit 
tion 

Uo = free-stream velocity 

u = velocity at any point in the boundary layer 

Ue = velocity in the boundary layer at n = e 

U = limit of « at the outer edge of the boundary 
layer 

AH = total head at a point in the boundary layer 
relative to the corresponding value for fully 
developed turbulent flow 

6 = thickness of Pohlhausen laminar boundary-layer 
profile 

@ 
0 = momentum tnickness ff (u/U){1 — (u/U)|dn 
0 

n = 

Ne = 

x = n(U/vs)'!2 

v = kinematic coefficient of viscosity 

= Uoc/v 

R. = 

a,m = coefficients in Eq. (1) 

Qa, = coefficients in Eq. (11) 

K,X = form parameter and shape parameter, respec- 


tively, in the Karman-Pohlhausen solution of 
the laminar boundary-layer flow 

F(n), G(n) = Pohlhausen functions for the laminar boundary- 
layer profile 


(1) INTRODUCTION 


) pragesnd in surface roughness, insofar as it may affect 
the boundary layer, has been confined in the past 
almost exclusively to the artificial kind produced by 
such disturbances as rivet heads. With the introduc- 
tion of the low-drag airfoil, however, and the more re- 
cent developments in the practical control of the lami- 
nar layer by suction, the problem of natural roughening 
due to insects has arisen. As long ago as 1945, Smith 
and Higton,' when investigating reduction of profile 
drag, drew attention to this source of premature transi- 
tion, and found that, on a low-drag airfoil, the rough- 
ness in question can produce a rise of from 50 to 100 per 
cent over the effective incidence range. Other research 
workers concerned with full-scale flight experiments on 
boundary-layer suction have also referred to the difli- 
culty.”,* So far, however, no papers appear to have 


been published on research specifically designed to ex- 
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amine the aerodynamic properties of natural rough- 
ness, or their relation to those established for various 
artificial forms of surface irregularity, although Keeble,’ 
Atkins,* and Johnson® have measured the chordwise 


spread of insect traces on a number of aircraft. 

The present work is concerned with this less familiar 

form of surface roughening, first, from the point of view 

| of its composition, and, second, in relation to its effect 

upon the stability of the laminar boundary layer. 

Quantitatively, the former aspect involves features such 

as the total streamwise extent of the roughness, the mag- 

nitude and distribution of excrescence height, and the 

intensity——that is to say, the proportion of unit area of 

_ the surface covered by deposits—while the latter is es- 

irface, sentially the problem of determining the critical condi- 
ie tions for forward movement of the transition. 

Unlike artificial disturbances, which may be ar- 

ransi- ranged in any conceivable pattern of a regular or irreg- 
ular character, those produced naturally have an in- 
variant distribution consisting of a confined area near 
the leading edge where the excrescence height reaches a 
ndary very large value, followed by a much more extensive 
region of considerably smaller elements. Normally, the 
Pina severe initial roughness is sufficient to bring the transi- 
ully 
tion right forward so that virtually the whole of the 
layer boundary layer becomes turbulent, but there is also 
another critical condition of practical significance asso- 
dn ciated with the limit at which the roughness ceases to 
constitute a disturbance to the passage of a laminar 
layer. This is investigated in Sections (6) and (8). 

It need hardly be observed that, in the kind of work 
under discussion, the problem of developing suitable 
methods of research has to be carefully considered. 
Clearly, the best course, from the experimental point of 
view, is to adopt a method which depends on direct ob- 

spec- servation in flight. Not only does this provide all the 
m of relevant full-scale conditions, but ensures a proper 
sampling of the aerial population responsible for 
a: the roughening of the surface. At the same time, it is 
equally obvious that, if the many parameters involved 
are to be adequately investigated, this procedure will be 
prohibitively slow and costly, particularly when it is 
fect recalled that, in most localities, favorable meteorological 
vast conditions only occur during a comparatively short pe- 
by | riod of the year. Moreover, the method lacks inher- 
luc- | ently the degree of control possible in the laboratory. 
re- Alternatively, suitable wind-tunnel experiments have 
mi- been devised for the investigation of the roughness pro- 
ing file [see Sections (2) and (3)] and the boundary-layer 
ith stability [see Sections (5) and (6)]. On the other hand, 
file for the reasons outlined in Section (3), these experi- 
nsi- ments do not strictly indicate the total streamwise ex- 
gh- tent of the roughness appropriate to the actual flight 
per’ conditions represented in the tunnel. Here, however, 
rch | theory is of considerable value. Thus, on the assump- 
on tion that an insect reacts like an inorganic particle to 
fi- the velocity field set up by the passage of a body (an 
ive assumption for which reasonable arguments may be ad- 
eX- vanced), it is possible to solve the hodograph equa- 
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tions, given the resistance law for the particle.* Hence, 
particle paths and impact velocity—.e., particle veloc- 
ity normal to the surface at a point of contact-—may be 
calculated. The limit of roughness then follows from 
the fact that this must occur where the local impact 
velocity becomes equal to the rupture velocityt of the 
most fragile particle encountered—a condition which 
may be assumed to arise very frequently at the low al- 
titudes associated with take-off and landing. Esti- 
mates of the total extent of roughness based on some of 
the more comprehensive flight experiments recorded in 
reference 5 (virtually the only reliable and adequate 
data which have been published) are in good accord 
with the measurements. 


As a result of this work, it has also been possible, on 
certain tentative, but plausible, assumptions regarding 
the impact conditions, to predict the maximum ex- 
crescence height which will normally occur at any 
point on a fully roughened surface. Agreement with 
the present observations [see Section (3) | is very satis- 
factory outside the confined and heavily contaminated 
region near the leading edge, but, as shown in Section 
(8), this restriction is of no consequence with respect to 
the conditions of practical interest investigated in Sec- 
tion (6)—namely, when the initial roughness, or at 
least its effect upon the boundary layer, is eliminated. 
In this case, therefore, the above analysis, in combina- 
tion with the present data relating to the criteria for 
stability of the boundary layer, may be used to predict 
the extent of the roughness zone in which a laminar 
layer cannot be self-sustaining. Again, calculation 
checks well with the experimental data. Conse- 
quently, although it has been thought preferable to 
present these various comparisons with the details of 
the theoretical investigation in the separate paper re- 
ferred to previously, it may be said now that the con- 
sistency between the two lines of approach affords good 
evidence for the validity of the general conclusions. 


(2) EXPERIMENTAL TECHNIQUE FOR THE GENERATION 
OF NATURAL ROUGHNESS IN THE LABORATORY 


To investigate the degree and distribution of rough- 
ness, as expressed by the excrescence height ¢, and to ex- 
amine its effect upon the behavior of the boundary 
layer, resort may be made to simple wind-tunnel ex- 


periments. So far as the generation of the roughness is 


* Fortunately, from Professor Hocking’s experiments,® the 
drags of insects over a wide range of size, as measured by a deli- 
cate torsion balance in a smail wind tunnel, are known. 

7 Rupture velocity has been defined in the present research as 
the impact velocity at which an excrescence not exceeding 3 
microns thick is formed. This was about the smallest displace- 
ment which could be measured consistently [see Section (3)]. 
In the case of flies from six different families, including the fruit 
fly, Drosophila, and the common housefly, which have been used 
in the present investigations, wind-tunnel measurements have 
shown extreme deviations from the mean values of about +2 
ft. per sec. Field experiments to determine the minimum rup- 
ture velocity gave a mean value with possible variations of +4 


ft. per sec. 


git 
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Fre. 1. 


Wind-tunnel installation for the roughening of an 
airfoil. 


concerned, it is obvious that the process must follow 
that which occurs in flight, and consequently, some pro- 
cedure by which live insects may be discharged upon the 
surface is necessary. This in turn raises the question of 
how the very numerous and diverse aerial population 
may be adequately represented in the laboratory, but 
fortunately both theory and experiment indicate that, 
over the region where « becomes negligible insofar as it 
can cause any further disturbance to the boundary 
layer,* the kind of insect, within the relatively fragile 
group essentially responsible for roughness, is unim- 
portant. Practically all the work considered in the 
present paper, therefore, has been performed with the 
fruit fly, Drosophila, since it is very simple to breed in 
abundance, and, at the same time, has a rupture veloc- 
ity of about 46 ft. per sec., which is not greatly differ- 
ent from the minimum measured in the field— namely, 
35 ft. per sec., approximately. It is, therefore, suitably 
representative, though in fact, it forms a relatively small 
group of the aerial population. By contrast, a few ex- 
periments have been made with the housefly, which is 
rather more than twice the size of the fruit fly, and has 
a rupture velocity of 39 ft. per sec. 

The simple device used for discharging insects into 
the airstream consisted of a Perspex tube 6 in. long, 
with an outside diameter of 1 in. and a bore of 1/2 in. 
At each end, a brass disc was fitted. These discs were 
soldered eccentrically to a common spindle in the wall 
of the tube, so that, by rotating the discs, the openings 
were sealed, or exposed, simultaneously. A stand con- 
sisting of an adjustable vertical pillar and horizontal 
plate, to which the tube could be clamped in any de- 
sired position, completed the instrument. 

To roughen a surface, the tube was mounted in the 
tunnel, with its axis into wind, at some suitable dis- 
tance upstream of the model. It was then charged with 
a number of live insects (commonly between 50 and 
100) and sealed. Finally, when steady flow conditions 
in the tunnel were established, the insects were dis- 


* In other words, with respect to the particular flow conditions 
considered in Section (6). 


charged by rapid opening of the tube, the latter opera- 
tion being performed with the help of a cord running 
from a small lever on the upstream disc to a point out- 
side the tunnel. Fig. 1 shows the general arrangement 
for the roughening of an airfoil. In this manner, any 
desired extent of surface could be gradually treated by 
successive displacements of the tube. For every set- 
ting, however, a number of discharges were normally 
required before the local roughness was fully estab- 
lished. 

A typical example of roughness produced in a wind 
tunnel by fruit flies is presented in Fig. 2. It refers toa 
13.5 per cent thick symmetrical airfoil at an incidence of 
7°, and was recorded on stout paper stretched tightly 
over the surface of the section. For comparison, Fig. 3 
is a photograph of natural roughness obtained by John- 
son on the Armstrong-Whitworth A.W.52 from his ob- 
servations in flight.® 


(3) MEASUREMENT OF EXCRESCENCE HEIGHT 


The magnitude and distribution of excrescence height 
e were studied on a plane surface and on two airfoils— 
namely, NACA 66-009 and NACA 66;-018 (see refer- 
ence 7). The plate was made of 3/4-in. reinforced ply- 
wood with a duralumin covering, and the unroughened 
surface tapered, to form a sharp leading edge, over a 
distance of 6in.fj Each airfoil consisted of a sheet dur- 
alumin profile supported on a structure of wooden spars 
and ribs, the chord, as in the case of the flat plate, being 
5 ft. 

7 For the reason given in the last paragraph of Section (3) 
below, only the lower surfaces of the plate and the airfoils were 
roughened. 
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Fic. 2. Example of roughness produced in the wind tunnel with 


fruit flies (Drosophila) on the surface of an airfoil of thickness- 
chord ratio 0.135 at an incidence of 7°. 
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CHARACTERISTICS 


Experiments were conducted in the 7 X 5 ft. low- 
speed tunnel® at Brough with the sections mounted be- 
tween the floor and roof on incidence pivots, and with 
the discharge tube placed approximately one chord up- 
stream of the leading edge (see Fig. 1). 

To obtain a record of the roughness which could be 
conveniently examined with a microscope, a soft cop- 
per strip, 1/4-in. wide, was stretched over the surface at 
the median section of the airfoil and secured with cello- 
phane tape along each edge, as well as at both ends, to 
ensure that it followed the profile accurately. Thus, 
after the strip had been roughened, as described in Sec- 
tion (2), it was removed and mounted on a specially 
prepared rule for easy manipulation under the micro- 
scope. Measurements of « were then made at inter- 
vals along the strip corresponding to 0.01 chord ap- 
proximately, and at each station several transverse 
readings were taken, varying in number from 1 to 12, 
according to the extent of the local roughness. The 
procedure to obtain ¢€ consisted, first, in focusing the 
microscope on the minute scratches in the surface of the 
copper, and secondly, on the limit of the roughness at 
the center of the field of view (0.012 in. dia.). Finally, 
a check of the surface reading was taken. Under the 
maximum magnification—namely, 400—repeated meas- 
urements of « varied by 2 to 3 microns. Further, be- 
fore the main experiments were started, several strips 
were roughened on the NACA 66-009 airfoil at the same 
wind speed and angle of incidence to determine whether 
the significant feature of the roughness, namely, the 
envelope to ¢ [see Section (4)], recurred consistently. 
This was found to be essentially the case. 

It is not intended to discuss here the influence of air 
speed on the formation of the roughness, since this fea- 
ture is more suitably dealt with in the theoretical in- 
vestigation referred to in Section (1). It will be suf- 
ficient, therefore, to say that the chordwise width and 
shape of the excrescence envelope just mentioned are 
both affected. Again, as we have already observed in 
the Section above, the extent of the roughness produced 
in a wind tunnel is not precisely the same as that which 
occurs in flight on an identical surface at the same inci- 
dence and wind speed. This is due to the fact that the 
relative velocity between particle and aircraft in steady, 
level flight remains constant if we make the legitimate 
assumption that the small random motion (free activ- 
ity) of the former may be neglected, and ignore, for 
the moment, the effect of the induced velocity field. 
On the other hand, in the laboratory experiment the 
particle is accelerated from rest by an air stream of 
constant velocity. Thus, unless it is discharged 
from a point far upstream of the airfoil, it will strike the 
surface in a state of acceleration and at a speed less than 
that corresponding to the limiting condition of zero 
slip, the residual acceleration and deficiency of speed 
being clearly a function of the distance between the dis- 
charge tube and the airfoil. It is also obvious that su- 


perposition of the induced velocity does not modify 
this difference between the flight and laboratory condi- 
tions in any fundamental way. Hence, for any wind 


OF ROU 


INSECTS 267 


FROM 


GHNESS 


Fic. 3. Natural roughness generated in flight on the wing of 
the Armstrong-Whitworth A.W.52, as recorded by D. Johnson.® 


tunnel of normal proportions, the impact velocity, and 
therefore the extent of the roughness, will, at a given 
wind speed, be less than their respective values at the 
same relative speed in level flight. But although sam- 
ple estimates have been made to obtain the order of the 
discrepancy, it requires an extensive exploration of the 
flow, together with prohibitively lengthy calculations, 
to establish the correction for even one combination of 
the parameters involved. Fortunately, in the immedi- 
ate investigation we are concerned exclusively with the 
critical value of ¢ in relation to chordwise position and 
the associated pressure distribution, to which end it is 
not essential that the roughness patterns shall accord 
strictly with those under flight conditions. In other 
words, the necessary data for the purpose of a gen- 
eral analysis [see Section (8)] may be adequately ob- 
tained from arbitrary profiles produced, for example, 
by varying the incidence at a constant wind speed,* 
the procedure adopted in the present experiments. A 
fairly high speed, corresponding to a chord Reynolds 
Number of approximately 7 X 10°, was chosen to make 
the roughness as extensive as possible. Even so, the 
difficulty of meeting this requirement adequately on 
the upper surface becomes acute. It is true that at 
zero incidence both surfaces are roughened by equal, or 
very nearly equal, amounts, depending on whether the 
airfoil is symmetrical or slightly cambered, but as the 
incidence is increased, roughness on the upper surface 
diminishes rapidly to a vanishingly narrow strip adja- 
cent to the leading edge, thus making exploration of the 
boundary-layer flow, as described in Sections (5) and 
(6) virtually impossible to perform with any precision. 
In the case of the lower surface, on the other hand, the 
opposite condition holds, and the above difficulty does 
not arise. It appears inevitable, therefore, that ob- 
servation must be confined to the latter region, but it 
seems unlikely that this restriction will be of any con- 
sequence in a practical problem to which the data, as 

*It need hardly be said that this must exceed the minimum 
wind speed to give rupture. 
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analyzed and checked in Section (8), may be applied. 
Accordingly, all results which we shall now discuss re- 
fer to conditions on the lower surface. 


(4) MAGNITUDE AND DISTRIBUTION OF EXCRESCENCE 
HEIGHT 


Figs. 4 and 5 show typical distributions of « in the 
case of NACA 66;-018. The measurements of local ex- 
crescence height are represented as vertical lines sur- 
mounted by a horizontal mark. The profile obtained 
by drawing a curve through the extremities of the or- 
dinates is, without exception, exceedingly irregular, and 
differs in detail from one specimen to another under the 
same test conditions. Nevertheless, despite the ran- 
dom character of the individual readings, there is evi- 
dence of two fairly consistent features. First, the de- 
posits tend to accumulate locally in groups which 
may be regarded as larger, compound excrescences 
with moderately well-defined contours; second, these 
accretions occur at intervals in the streamwise direc- 
tion and diminish regularly in height to zero at the ab- 
solute limit of the roughness. Thus, it would seem rea- 
sonable to infer that the surface flow is principally dis- 
turbed by these major irregularities, in which case the 
distribution of significant excrescence height will be 
given by the curve drawn through the summit of each 
mass. This hypothesis appears to be well supported 
by the results of the analysis in Section (8). 

From the envelopes so obtained, other features of 
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the roughness are immediately apparent. Particularly 
noticeable is the very rapid increase of excrescence 
height in a concentrated area adjacent to the leading 
edge, the maximum value of ¢« being located in the vi- 
cinity of the forward stagnation point. The same con- 
clusion was also reached by Johnson’ in his flight ex- 
periments. Downstream of this region, however, ¢ 
is much diminished, and, within broad limits, the en- 
velope decreases linearly to zero. 

Fig. 6, in conjunction with Figs. 4 and 5, shows the 
effect of incidence on the distribution of « (only the en- 
velopes are presented for clarity), while Figs. 7 and § 
illustrate, respectively, the influence of an appreciable 
change in the particle characteristics, and of surface 
curvature and airfoil thickness. The last two parame- 
ters clearly have a substantial bearing on the genera- 
tion of the roughness. Their separate effects are not 
easily resolved experimentally, but broadly speaking, a 
thickening of the section, and hence an increase of the 
curvature, reduces the total streamwise extent of ‘the 


roughness, the amount of the reduction being much | 


greater for thicknesses between 9 and 18 per cent of the 
chord than between 0 and 9 per cent. In contrast, a 
marked reduction of ¢, especially near the leading edge, 
takes place when the airfoil tends to the flat plate pro- 
file. These conclusions are fully confirmed by theory, 
and, as the examples I, II, and III in Section (8) dem- 
onstrate, a profile having a thickness of, say, 3 to 4 per 
cent of the chord is, even at the largest full-scale Rey- 
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Fic. 4 (Jeft). Distribution of roughness obtained in the wind tunnel with fruit flies on the lower surface of NACA 663-018; high-speed 
flight condition, R. = 6.9 X 16°. Fic. 5 (right). Corresponding roughness distribution on the lower surface of NACA 663-018 under 
take-off or climb conditions approximately; R, = 6.9 X 10°. 
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Fic. 6 (left). Effect of incidence on the lower surface roughness envelope; NACA 66-009, R. = 6.9 XK 10% Fic. 7 (right). Effect of 
a pronounced change in the particle properties on the lower surface roughness envelope; NACA 66-009, R. = 6.9 X 10°. 
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Fic. 8. Effect of surface curvature (wing profile) on the lower surface roughness envelope; R, = 6.9 X 10°, 


nolds numbers, much less susceptible to natural rough- 
ening than one of more usual proportions. 

Change of incidence has a very obvious effect. As 
Figs. 4, 5, 6, and 16 clearly show, the total extent of the 
roughness, x,, on the lower surface varies almost lin- 
early with incidence (a conclusion which is supported 
theoretically), and increases as incidence increases. The 
converse holds on the upper surface, though in this case 
the variation predicted by calculation is not linear. 

The way in which the properties of the particle—.e., 
its size, mass, drag, and rupture velocity—influence the 
distribution of the roughness is not immediately obvi- 
ous. Analytically, it is found to depend on the nondi- 
mensional ratio 


= (1/2)Cp(p/pi) (¢/2) 


where Cp is the drag coefficient of the particle, p the 
density of the surrounding fluid, p; the density of the 
particle, c any typical length of the body, and / a cor- 
responding length of the particle. In the case of a wing, 
when c is taken as the local chord, and the particle is an 
insect of body length /, the possible limits of « for any 


normal size of aircraft are approximately 0.7 < « < 3.5, 
Over this range, calculation shows that, on the lower 
surface, the extent of the roughness and the rearward 
distribution of ¢ are only slightly affected. This conclu- 
sion is reasonably well supported by the results of Fig. 
7, the value of x in the experiments being 0.9 for the 
fruit fly and 0.4 for the housefly.* The pictorial evi- 
dence of Figs. 9(a) and 9(b), obtained from photo- 

* The low values of « in the above experiments compared with 
the values quoted for full-scale conditions are due to the small 
wing chord—namely, 5 ft. 


Fic. 9(a). Cross section of the roughness from fruit flies near 
the leading edge of an airfoil. 
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Fic. 9(b). cross section of. roughness from 
houseflies. 
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Fic. 10(a). China clay pattern behind an isolated excrescence 
(fruit fly) at approximately 4 per cent of the chord from the 
leading edge. 
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Fic. 10(b). Similar patterns behind a cluster of excrescences 
(fruit flies) located between the leading edge and 2.5 per cent of 
the chord. Note absence of turbulent wake in Pi case of the 
excrescence at the leading edge, and gradual development of such 


wakes over the initial 1 per cent of the chord. 


graphs taken more or less tangentially to the surface 
near the leading edge of an airfoil used in some early 
exploratory experiments, again confirms the later 
measurements, and shows clearly the nature of the 
deposits for the same two insects. In the case of the 
upper surface, on the other hand, theory predicts a 
moderate variation of the extent of the roughness with 
x, but it has not been possible to check this conclusion 
experimentally owing to the difficulty discussed in 
Section (3). Nevertheless, the nature of the particle, 
within the limitations relevant to the natural rough- 
ening of surfaces, does not appear to be an important 
factor so far as the distribution of ¢ is concerned. 


(5) PRELIMINARY INVESTIGATION OF THE FLOW OVER 

THE SMOOTH SURFACE, AND OF THE MANNER IN WHICH 

DEGREE OF ROUGHNESS AFFECTS CONDITIONS IN THE 
BOUNDARY LAYER 


The state of flow in the boundary layer, as affected by 
roughness, was examined on the lower surface of the 
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NACA 66-009 airfoil with a small hypodermic pitot 
tube flattened at the exposed end to a rectangular open- 
ing measuring 0.003 in. normal to the surface by 0.045 
in. wide. The length of the tube along wind was also 
about | in. 

The mechanism for the tube displacement normal to 
the surface consisted of the usual micrometer system 
with a screw of 0.05 in. pitch and a vernier wheel ca- 
pable of reading to0.001in. The latter was mounted on 
the outside of the tunnel in a position convenient for ad- 
justment of the micrometer during an experiment, the 
zero position on the surface being recorded electrically 
by contact of the lip of the pitot tube with a flush-fitted 
stud in the airfoil. Corrections to total head and dis- 
placement were taken directly from the calibration data 
for instrument A in Tables VII and VIII of reference 9, 
since the dimensions relating thereto were almost the 
same as those of the tube used in the present research. 

In addition, static pressures were measured at a num- 
ber of chordwise surface holes on the forward part of the 
airfoil (0 < x < 0.25c), the readings, as in the case of the 
boundary-layer surveys, being recorded on a wide-range 
Chattock manometer for angles of incidence of 0°, 1°, 
3°, and 6°. Finally, all measurements were made at a 
nominal wind speed of 250 ft. per sec., which, for the 
associated stream temperature, during a prolonged run, 
of approximately 50°C., corresponded to a chord Rey- 
nolds Number of rather less than 7 X 10°; i.e., the same 
as that in the investigation of excrescence height [see 
Section (3)]. The data about to be discussed also re- 
fer to observations at the median section of the airfoil 
where the flow was closely two-dimensional. 

Before proceeding to the investigation of the effects 
of roughness, it was necessary to establish the boundary- 
layer conditions on the bare, smooth surface. In or- 
der to give good contrast to the deposits, the former was 
sprayed with three coats of white cellulose paint and 
polished with Durex paper to an egg-shell finish. Bound- 
ary-layer surveys were then made for each of the above 
angles of incidence at several chordwise positions, 
namely, at 0.05c, 0.10c, 0.15c, and 0.20c. In every 
case the profile was characteristically laminar, and val- 
ues of the displacement and momentum thicknesses, as 
determined from the experimental data, were in excel- 
lent accord with those calculated from the static pres- 
sure distributions by the method of ref. (10). Addi- 
tional experiments to detect the position of transition 
by the china clay technique (11) indicated that this oc- 
curred at between 30 and 40 per cent of the chord from 
the leading edge within the incidence range investi- 
gated. In both the case of the airfoil at present under 
discussion (NACA 66-009), and of NACA 66;-018, with 
which we shall be concerned in Section (6), it is clear, 
therefore, from Fig. 16, that the maximum extent of the 
roughness lay within the region of laminar flow on the 
smooth surface. 

Again, an essential preliminary to the main theme 
of the experimental investigation arose from the need to 
examine the way in which the mean intensity of rough- 
ness—i.e., the proportion of area covered by ex- 
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crescerices on a small finite element of the surface—in- 
fluences the behavior of the boundary layer. But 
apart from the necessity to appreciate quantitatively 
the significance of this feature for the correct interpreta- 
tion of data relating to arbitrary degrees of roughening, 
such information is not without direct interest. 

Initial china clay records of the wakes generated by 
individual fruit flies and groups of the same (Fig. 10) 
showed, as may be expected, features typical] of those 
observed with small pin heads and other shapes of 
pimple which have been used in work on artificial rough- 
ness. On completion of these qualitative experiments, 
the airfoil was roughened in stages, and at each of the 
incidences previously mentioned—that is, 0°, 1°, 3°, 
and 6°—a boundary-layer traverse was made at a point 
well downstream where the surface was known, from 
experiment, to be aerodynamically smooth.* Thus, by 
this gradual build-up of the roughness, it was possible to 
follow the changes in the boundary-layer flow. To meas- 
ure the intensity, a series of photographs of the af- 
fected surface was taken. These records resembled, of 
course, the example shown in Fig. 2. From prints, 
considerably enlarged, the dark and clear areas, indi- 
cating, respectively, excrescences and smooth surface, 
were then cut out jig-saw fashion and weighed in 
separate groups on a chemical balance. Hence, as- 
suming uniformity of weight of the photographic paper, 
the area of surface which had been roughened could be 
determined. This very laborious, and admittedly ap- 
proximate, procedure was only adopted in the heavily 
contaminated area near the leading edge. Further 
downstream, where there is a considerable thinning 
out, and a fair uniformity of size and shape of the ex- 
crescences, the roughness areas were calculated from 
readings with a measuring microscope. Intensities 
were also determined for intervals of 1 per cent of the 
chord over a distance of 4 in. transverse to the flow, 
this being broadly the width of uniform roughening ob- 
tainable at one spanwise position of the discharge tube. 
All experiments were performed with the fruit fly as 
roughening agent. 

Fig. 11 is a typical example of the way in which the 
boundary layer changes from the laminar state to that 
of fully developed turbulence as the roughness is gradu- 
ally increased in intensity. It will be observed that 
the transitional process takes place correspondingly at 
a fairly steady rate, and is very readily promoted. 
Thus, for the conditions of Fig. 11, a perceptible change 
of profile toward the turbulent form has occurred when 
the roughness has spread to no more than 2 per cent of 
the chord, with an intensity near the leading edge of 
only 5 per cent. Moreover, it appears to be virtually 
complete when the roughness has increased to about 65 
per cent of its maximum chordwise extent, and the 


* The downstream limits of the roughness allowed these sur- 
veys to be made at the chordwise positions chosen in the experi- 
ments on the smooth surface—namely, at 0.05c, 0.10c, 0.15c, and 
0.20c—corresponding, respectively, to angles of incidence of 
0°, 1°, 3°, and 6°. It thus became unnecessary to remeasure the 
velocity distribution for the unroughened case. 


leading-edge intensity to approximately 40 per cent. 
Any further increase of intensity had no effect on the 
velocity profile which continued to be that designated 
“fully roughened surface” in Fig. 11. It accords very 
well with the logarithmic law of turbulent boundary- 
layer theory, and there appears to be no doubt, there- 
fore, that the flow changes from the purely laminar to 
the completely turbulent state during the indicated 
growth of the roughness. Similar results were ob- 
tained at the other angles of incidence. It would seem, 
consequently, that, in flight under conditions favorable 
to natural roughening, sufficient insects may be trapped, 
as instanced in Fig. 3, to destroy completely an other- 
wise laminar flow over the forward part of the wing. 


(6) EXPERIMENTAL DETERMINATION OF CRITICAL 
EXCRESCENCE HEIGHT 


As shown in Figs. 4 to §, a recurrent feature of natu- 
ral roughness is the very pronounced peak in the en- 
velope of excrescence height near the leading edge. 
The data suggest that this intense source of disturbance 
will be sufficient to bring transition forward so that the 
boundary layer becomes turbulent over virtually the 
whole of the roughened surface. In this Section, ex- 
periments are described to investigate this point. As 
we shall see, the results show that transition does, in 
fact, occur in the region concerned, even at the lowest 
Reynolds Number considered—namely, at a value, 
based on the chord, of approximately 3.5 X 10°. 
Consequently, since the local magnitude of the ex- 
crescence is sensibly constant, being composed mainly 
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Fic. 11. Example of the alae os roughness intensity on the 
boundary-layer flow. Survey rey" on the lower surface of 
NACA 66-009 at x = 0.05c, a = 0°, and R. = 6.9 X 105. 
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Fic. 12. Example of the behavior of the boundary layer as 

the roughness is gradually removed in the downstream direction. 

Survey made on the lower surface of NACA 66-009 at x = 
0.05c, a = 0°, and R. = 6.9 X 10. 


of more or less whole insects, it is evident that this col- 
lapse of the laminar flow will hold with even greater 
force in the range of full-scale Reynolds Numbers, and, 
as a result, any attempt to reduce the drag of the smooth 
surface by means of boundary-layer control will be 
very largely, if not completely, nullified. For this rea- 
son, various tentative methods of eliminating the sig- 
nificant roughness have been proposed—e.g., by a tem- 
porary and jettisonable cover around the leading edge 
of the wing, or by a mechanical scraper which removes 
the excrescences—but the usefulness of such devices 
may be shown to depend very much on the amount of 
surface which has to be protected or cleaned. Conse- 
quently, it becomes of some practical interest to deter- 
mine the downstream limit at which the surface, though 
still rough, may again be regarded as aerodynamically 
smooth; or in other words, to establish the critical 
value of ¢ defining this section. 

This aspect of the investigation was studied by suc- 
cessively removing the roughness in a series of narrow 
spanwise strips. Thus, starting at the leading edge, a 
small chordwise width of the original smooth surface 
was exposed, thereby restoring laminar flow across the 
element. Clearly, however, transition recurs at the 
new upstream boundary line of the roughness if ¢« at 
this position is equal to, or greater than, the local criti- 
cal value.* Hence, by gradually cutting back the 
roughness to uncover more and more of the smooth sur- 
face, a point is eventually reached at which the bound- 
ary layer continues to remain laminar. This, then, 
gives the chordwise location and critical value of ¢€ ap- 


* Unlike the two-dimensional disturbance, such as a transverse 
wire, which brings the transition forward gradually, the three- 
dimensional equivalent, of which natural roughness is an ex- 
ample, induces transition immediately at the critical excrescence. 
The mechanism of transition in the two cases is discussed further 
in Section (7). 


propriate to the particular conditions of flow over a 
roughened surface at present under discussion. 

Observations were made on the lower surface of both 
the NACA 66-009 and NACA 66;-018 airfoils. In the 
case of the first, the roughness was mainly removed in 
spanwise strips having a projected width along the 
chord of 0.3 in., though, where circumstances per- 
mitted, the interval was increased to double this value. 
The effects of these successive eliminations upon the 
boundary layer were then recorded from pitot surveys 
at the same chordwise positions and corresponding 
angles of incidence as in the earlier work on roughness 
intensity [see Section (5) }. 

Fig. 12 is a record of the results at zero incidence. 
The gradual change from the fully developed turbulent 
profile associated with total roughness to the laminar 
distribution when the roughness over approximately 2 
per cent of the chord from the leading edge has been re- 
moved is very clear. Similar results were obtained for 
the other angles of incidence. Thus, the point at which 
the excrescence height becomes critical is indicated im- 
mediately from each family of profiles. To estimate 
its value precisely, the momentum thickness 6 for every 
individual curve in a group was determined, and the 
results were then plotted as a function of the extent Ar 
of the roughness removed—.e., of the aggregate of all 
strip widths up to the point concerned. The relevant 
curves are shown in Fig. 13, from which it is obvious 
that the distance x, at which the roughness has become 
insignificant occurs at the point where @ has fallen to the 
value for the laminar boundary layer. 

Although the above procedure clearly provides an 
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Fic. 13. Determination of x; from plots of 6, as given by the 
measured boundary-layer profiles, for the lower surface of NACA 
66-009; R. = 6.9 X 10° 
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Fic. 14. Details of the fixed pitot probe used for the determina- 
tion of transition on NACA 663-018. 


instructive and detailed account of the behavior of the 
boundary layer, it is obviously very laborious. Conse- 
quently, having obtained so complete a record in the 
case of NACA 66-009, it seemed scarcely worth while 
to repeat the process on NACA 66;-018. Instead, the 
more rapid technique of the fixed pitot probe, or rather 
of the double arrangement described in reference 12, 
was employed for the location of transition. Fig. 14 
gives details of the instrument and of its location rela- 
tive to the surface. The pitot tubes proper were 
formed from quartz glass capillaries cemented with Ar- 
aldite into steel tubes ((0.104-in. external diameter) 
which were soldered to brass angle brackets, the latter 
being aligned with the wind direction and secured to the 
wing surface asshown. Thus, measuring total head, for 
convenience, relative to the local value in the fully de- 
veloped turbulent boundary layer, it is clear, from com- 
parison of the laminar and turbulent profiles (see Figs. 
11 and 12) with the pitot tube settings given in Fig. 14, 
that, when transition takes place, the inner tube will 
record a negative pressure increment and the outer a 
positive increment. 

Fig. 15 shows some of the transition records for 
NACA 663-018. Only the extreme angles of inci- 
dence—namely, 0° and 6°—were considered in the 
case of this airfoil, but observations were made at a 
chord Reynolds Number of about 3.5 X 10° as well as 
7 X 10% Apart, however, from the use of the fixed pi- 
tot probe instead of the pitot tube and micrometer, the 
experimental procedure was precisely the same as in the 
case of NACA 66-009—that is, strips of the roughness 
0.3- to 0.6-in. wide in the direction of the chord were re- 
moved successively until the critical value of € was 
reached. Also, the probe position at each angle of in- 
cidence was the same as for the traversing pitot in the 
previous experiments. Actually, rather more of the 
roughness was removed than the amount upstream of 


the critical point in order to make quite sure that the 
laminar profile had been restored. This is indicated in 
Fig. 15 by the constancy of AH after the transition 
jump. The distance x, at which e became critical was 
then judged to occur at the value indicated by the 
broken vertical line on each graph. Finally, when x, 
has been established by either of the techniques dis- 
cussed above, €,,;: follows immediately from the data of 
Figs. 4, 5, and 6. 

Curves of the velocity at the outer edge of the bound- 
ary layer (U’), the total chordwise extent of the rough- 
ness (x,), the corresponding extent of the significant 
roughness (x,), and the critical value of €, as measured 
on the lower surfaces of the NACA 66-009 and 663-018 
sections, are presented in Fig. 16. For the reason 
given in Section (3), wind-tunnel measurements of x,, 
based on the technique of Section (2), do not yield re- 
sults which, apart from the possible effects of wall con- 
straint, are strictly valid for the same wing chord, speed, 
and incidence under flight conditions. However, pro- 
vided the discharge tube is placed at a sufficient dis- 
tance upstream of the leading edge to ensure that all 
particles are moving with substantially uniform and 
equal velocity as they enter the induced field of flow set 
up by the airfoil—a condition moderately well satisfied 
in the experiments—the discrepancy will be small. Con- 
sequently, comparisons of x, and x, are not seriously af- 
fected, and it is immediately evident from Fig. 16 that, 
for the relevant range of Reynolds Number at any rate, 
a large amount of the roughened surface causes no dis- 
turbance to a laminar layer passing over it. This point 
is extended to full-scale conditions in Section (8), where 
it is shown that the above conclusion holds up to the 
highest Reynolds Numbers attained in flight. 


(7) A MetrHop oF ANALYZING THE EXPERIMENTAL 
DaTA 


It is reasonable to suppose that the action of natural 
roughness upon the laminar boundary layer is similar to 
that observed with artificial disturbances of a three-di- 
mensional character, in particular of the sandpaper vari- 
ety, and consequently, it is of interest to examine the 
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Fic. 15. Determination of x; on the lower surface of NACA 
663-018 by means of the fixed pitot probe; R. = 6.9 X 10°. 
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results of the investigation described in the preceding 
Section from this point of view. 

Analyses of the experiments performed with indivi- 
dual surface projections—e.g., discs and cones—and of 
the same projections in a variety of patterns, for differ- 
ent pressure gradients, are given in references 13 to 20, 
those of references 16 and 20 being especially thor- 
ough. Accordingly, we shall consider here only those 
aspects of the research concerned which have a direct 
bearing on the present work. In the case of three-di- 
mensional elements located upstream of the point 
where transition occurs naturally on the aerodynami- 
cally smooth surface, it has been observed’ ?! that no 
significant effect upon the laminar boundary layer oc- 
curs until a critical point is reached at which transition 
moves forward catastrophically to the disturbance. 
This condition is associated with the appearance of tur- 
bulent spots which are generated by the projections 
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and discharged downstream. Moreover, it is found to 
be connected in a fairly well defined manner with the 
roughness Reynolds Number R, = «u,/v (where ¢ is the 
excrescence height, “, the velocity at a distance ¢e from 
the surface, and v the kinematic coefficient of viscosity), 
a criterion which was first advanced by Schiller.}?? 

The action on the boundary layer of three-dimen- 
sional particles is also in marked contrast to that of the 
two-dimensional irregularity, such as a wire transverse 
to the flow. In this case, it has been shown” that, if 5, 
is the distance of transition from the forward stagnation 
point on the smooth surface, and s, the corresponding 
distance in the presence of the ridge, then s,/so decreases 
progressively with increase of ¢/6,*, where 6,* is the 
displacement thickness at the projection, in a manner 
which is uniquely defined. Thus, the transition with 


¢ Schiller used the friction velocity uz = (70/p)!/? in place 
of u,., where 70 is the surface shear stress and p is the mass density, 
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Experimental values of U/U,, x;, xz and eérit on the naturally roughened lower surfaces of NACA 66-009 (top row) and 
NACA 66;-018 (bottom row). 
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this type of roughness gradually moves upstream from 
the position for normal laminar flow to a limit at the 
excrescence as the disturbance becomes more pro- 
nounced. 

With regard to the roughness produced by insects, 
the effect upon the boundary layer may be expected to 
show, as we have observed, a marked similarity to 
that of inorganic particles, particularly when the lat- 
ter are distributed at random over the surface. We 
shall examine, therefore, the results of the present in- 
vestigations on transition in terms of the roughness 
Reynolds Number 

For practically all cases,} the velocity u, at the peak of 
the excrescence lies within the laminar layer, and 
consequently, the relevant velocity profiles have to be 
established either experimentally or by calculation. 
Lack of the essential observational data has made it 
necessary to resort to theory in the analyses of artificial 
roughness, and the same procedure will be followed here, 
for although, as stated in Section (5), surveys of the lam- 
inar boundary layer were taken when checking the flow 
on the smooth surface, nevertheless, they were con- 
fined to the NACA 66-009 section, and to relatively 
large chordwise intervals. Consequently, the data are 
insufficient for a complete examination of the critical 
conditions at transition. At the same time, the avail- 
able measurements afford valuable checks of the calcu- 
lations (see Fig. 18). 

Two solutions of the laminar boundary-layer flow 
were considered and compared with experiment— 
namely, the Karman-Pohlhausen solution, as modified 
by Walz, and sunimarized in reference 18, and that of 
Falkner and Skan.** The latter is restricted to a par- 
ticular form of the streamwise velocity distribution at 
the outer edge of the boundary layer, but the formula is 
not so limited. On the other hand, if the velocity can- 
not be expressed analytically, it involves graphical dif- 
ferentiation. However, the roughness problem is es- 
sentially confined to small values of lift coefficient, and 
furthermore, is associated very largely with the low- 
drag type of airfoil. In which case, a satisfactory ap- 
proximation (see, for example, Fig. 16) proves to be the 
distribution studied by Falkner—namely, 


= a(s/c)”™ (1) 


where L’ is the velocity at the edge of the boundary 
layer, Uy is the free-stream velocity, s the distance along 
the surface from the forward stagnation point, c the 
chord, and a, m are constants. Eq. (1) may be directly 
applied to the Karman-Pohlhausen method to yield an 
appreciable simplification. Thus, substituting Eq. (1) 
in the integral relation for the momentum thickness 6, 


R.(6/c)? = 0.47(U/ d(s/c) (2) 


where R, is the Reynolds Number Uc/v, we obtain 


+ Exceptions may occur very close to the leading edge where 
the critical values of « can exceed the local thickness of the bound- 


ary layer. 


R(0/c)? = [0.47/a(5m + 1)] (s/c)'~" = ¢(a,m,s/c) 
(3) 


let us say. Hence, from Eq. (3) and the first derivative 
of Eq. (1), the form parameter A, which is defined by 
the expression, 


K = R(0/c)*d(U/ Uo) /d(s/c) 
reduces to 
K = 0.47m/(5m + 1) (4) 
But A is also related to the shape parameter 
= R(6/c)*d(U/ Uo) /d(s/c) 
by the relation 
K = [(37/315) — (4/945) — (A2/9072)|2\ (5) 
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Fic. 17. asa function of m for the Karman-Pohlhausen solu- 
tion of the laminar boundary-layer flow when U/U, is given 
by Eq. (1). 
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Fic. 18. Comparison of Karman-Pohlhausen laminar velocity 

profiles with measurements on NACA 66-009 when (a) a@ 

0°, s/c = 0.05, and (b) a = 6°, s/c = 0.10, at a value of R, 
6.92 X 10°. 
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Fic. 19. Analysis of the transition conditions for artificial 
roughness of a three-dimensional character. 


and so equating Eqs. (4) and (5), 
M(A) = 0.47m/(5m + 1) (6) 


where, for brevity, /(A) denotes the square of the quad- 
ratic in Eq. (5). 

Thus, for the velocity distribution defined by Eq. 
(1), A is purely a function of m, and is shown plotted in 
Fig. 17 over the interval of immediate interest. 

Again, Pohlhausen’s solution yields 


(0/6)* = f (7) 


where 6 is the boundary-layer thickness, and therefore, 
from Eqs. (3) and (7), 


JOURNAL OF THE AERO/SPACE SCIENCES—MAY, 1959 


Finally, 1, follows from the equation 

u./U = F(n.) +  G(n) (9) 
where F(n), G(n) have the well-known forms 

F(n) = 29 — + 

G(n) (1/6) (1 — 


Thus, for any observed value of ¢, the corresponding 
value of R, is immediately calculable from Eqs. (1), (3), 
(6), (8), (9), and the identity 


R, = (€/c) (U/Uo) (u/U) R. (10) 


the constants a, m in Eq. (1) being estimated to give the 
best fit to the experimental velocity distribution. The 
degree of approximation which this entails will be ap- 
parent from the typical examples shown in Fig. 16. 
Furthermore, it is found that deviations of a and m 
within the possible limits do not materially affect x,. 

The determination of the laminar velocity profile 
by the method of Falkner and Skan follows immedi- 
ately from Eq. (1), and is given by the series 


u/U = — > ka, x*™! (11) 


where, in the present notation, x represents the variable 
n(U/vs)'/*, and the coefficients a, are functions of m. 
The authors have evaluated a: by trial and error in the 
interval —0.09 < m < 2.0, and the remaining coeffi- 
cients can then be determined directly in terms of m 
and az. The relations for a; to a are to be found on p. 
7 of reference 23. Higher terms in the series (11) to 
give greater accuracy near the outer edge of the bound- 


ne = €/6 = (€/c) (¢/fR.) 7? (8) ary layer can be estimated, if necessary, for any particu- 
TABLE 1 
Fig. No. 
of Present a, 
Ref. Paper Symbol Surface Type of Roughness Re X 107% deg. Cx a m 
13 19 ok Airfoil 1, upper surface Spanwise row of circular, 5.0to8.8 — 0.347 1.420 0.0486 
cylindrical pegs spaced at 
3-in. intervals 
13 19 + Airfoil 2, upper surface 5.0to8.4  — 0 1.290 0.0245 
14 19 e 13 per cent low-drag Single conical pimple 4.8 0 = 1.350 0.1020 
section, upper surface 
14 19 A 5.8 0 — 1.350 0.1020 
20 19 © Flat plate, with 11 percent Single disc 6.5 to 7.6 1.119 0.0446 
trailing-edge flap —2 1.192 0.1004 
20 19 x Spanwise row of discs 5.3 to 7.6 1 1.094 0.0312 
spaced at 1.25-in. 0 _ 1.119 0.0446 
intervals 1.1388 0.0546 
20 19 A Flat plate, with 1l percent 1/4-in wide spanwise strips 5.0 to 7.9 1 - 1.094 0.0312 
trailing-edge flap of grit 0 - 1.119 0.0446 
—1 _ 1.188 0.0546 
—2 1.192 0.1004 
Present paper 21 © NACA 66-009 lower Insect (fruit fly) 6.92 0 - 1.158 0.0298 
surface 1 _ 1.193 0.0682 
3 _ 1.494 0.2092 
6 - 2.005 0.4370 
21 x NACA 66;-018 lower 3.46 0 - 1.454 0.1047 
surface 6 2.285 0.4320 
21 A _ 6.92 0 = 1.454 0.1047 
6 2.285 0.4320 
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lar value of m, but, apart from the comparatively rare 
instances indicated in the footnote on page 275 of the 
present paper, the first ten terms of the series normally 
proved sufficient for a satisfactory approximation. 

Values of u, were calculated by both methods in re- 
gard to all examples taken in the present survey of the 
data for artificial roughness, as discussed below. The 
Pohlhausen results proved to be consistently less than 
those from Falkner and Skan by a matter of 2 to 7 per 
cent. Further, from Fig. 18, which compares theory 
and experiment in the case of some typical profiles at 
approximately the transition points observed in the 
present investigation, it is evident that calculation 
slightly overestimates u in the middle and outer regions 
of the boundary layer where the majority of the critical 
conditions occur. Hence, the Karman-Pohlhausen 
procedure appears, on balance, to be in rather closer 
agreement with observation, and accordingly, the 
curves of critical Reynolds Number presented in Sec- 
tion (8) for both artificial and natural roughness are 
based exclusively on values of u, obtained by this 
method. 


(8) COMPARISON OF THE TRANSITION CONDITIONS FOR 
NATURAL AND ARTIFICIAL ROUGHNESS 


In view of the numerous analyses, to which we have 
already referred, concerning artificial roughness of a 
three-dimensional character, it is not proposed to con- 
sider here more than an extract of the results most 
closely connected with the data of the present investi- 
gation. Thus, Fig. 19 is a plot of Rone for a number of 
shapes and arrangements of projections summarized 
in Table 1. 

As will be seen, it is restricted to values of the chord 
Reynolds Number 2, essentially within the range se- 
lected for the study of transition on surfaces roughened 
by insects, namely, when R, lies between approximately 
5 X 108 and 9 X Over this interval, the variation 
of R, is confined to the limits indicated by the curves 
ABC and DEF. The area so defined is due principally 
to a recognizable dependence of R, on R,, the points 
nearer DEF referring in some cases mainly to the lower 
values of R,, and vice versa with respect to ABC. At 
the same time, random scatter of the data, and dif- 
ferences in the shape and distribution of the irregulari- 
ties, are contributory factors. The present analysis 
has also shown, though the evidence is omitted, that 
the correlation of results from the various investiga- 
tions is much less satisfactory when X, is less than about 
5 X 10.6 On the other hand, Johnson’s observations in 
flight'* at a chord Reynolds Number of approximately 
25 & 10° suggest that the curve ABC continues to hold 
very well for values of R, appreciably beyond a figure 
of about 10’, and hence, that R tends to become in- 


€ crit 


dependent of R, as full-scale conditions are approached. 
The analysis of the relevant data is not included in Fig. 
19, since the measurements relate to an irregularly 
spaced group of excrescences (paint spots), with a ran- 
dom distribution of «, and were designed to determine 


whether each individual pimple did, or did not, cause 
transition. Consequently, the critical curve of € can- 
not be definitely established, but, as Fig. 20 shows, the 
curve of €,,;, corresponding to the curve ABC of Fig. 19 
lies satisfactorily between the stable and unstable do- 
mains except very near the leading edge. However, as 
Johnson points out, the measurement of ¢ in this region 
is subject to an appreciable error, so that the compari- 
son may be markedly unreliable upstream of, say, 0.04c. 

The data of Fig. 19 also indicate a fairly well defined 
variation of RK, with s, which, for the small angles of 


€ crit 
incidence associated with roughness, does not appear to 
depend, in any consistent manner, on the streamwise 
velocity gradient. For this reason, an exhaustive rec- 
ord of R, |, in relation to changes of the last-mentioned 
variable has not been included in Fig. 19. 

The results of the corresponding analysis in the case 
of a naturally roughened surface are presented in Fig. 
21 which also includes re-plots of the curves ABC and 
DEF from Fig. 19. Subject to the rather limited nature 
of the evidence, it would appear that the same simple 
transition criterion, in terms of R holds for the kind 


€ crits 
of roughness produced by insects, and further, that it 
has closely the functional form appropriate to artificial 
types (see Fig. 19). Thus, at Reynolds Numbers of 107 
and more, the results of Fig. 21 indicate that the extent 
of natural roughness which is significant [see Section 
(6)] may be determined, in any specific case, on the 
basis of curve ABC. Consequently, given Uy, c, v, and 
U asa function of s, the value of 6 and the distribution of 
u with respect to 7, when s = constant, may be calcu- 
lated according to the procedure of Section (7). Hence, 
by interpolation, €,,;;, can be determined to satisfy the 


known value of R Finally, by plotting the curves 
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Fic. 20. Comparison of the calculated curve of eri: with the 


flight observations by Johnson (see reference 14, Part II). For 
the calculated curve, the same assumption was made as in refer- 
ence 14—namely, that the Vampire wing section could be suffi- 
ciently well represented by the 13 per cent thick airfoil investi- 
gated by Gregory and Walker at an incidence of 0°. 
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Fic. 21. Analysis of the transition conditions for roughness 
produced by insects (fruit fly). 


of ¢ defining the roughness envelope and of €,,;:, we ob- 
tain from the intersection = €,,;;. 

Illustrations of this graphical method of obtaining 
the transition point will be found in Fig. 23. The exam- 
ples have also been chosen to indicate the probable im- 
portance of natural roughness under the kind of full- 
scale conditions to which it may most commonly ap- 
ply should extensive laminar flow become a feature of 
aircraft design. In conclusion, therefore, we will con- 
sider the relevant calculations a little further. The es- 
sential details for the three cases considered are given 
in Table 2. 

These conditions, it may be noted, are not unrepre- 
sentative of, respectively, a modern feeder-line trans- 
port, a medium transport, and about the largest craft of 
this type when each is operating at its cruising speed. 
Thus, the relevant curves of €,,;; can be found immedi- 
ately from the data of Table 2, in combination with 
the method described above, for any specified distribu- 
tion of l’. It remains to establish the envelope of ¢ 
due to the roughening of the surface, or more particu- 
larly, of that part where the distribution is essentially 
linear, since the excrescence height adjacent to the 
leading edge is invariably significant [see Section (6) ]. 
In the former region, then, the equation for the enve- 
lope is, obviously, 


e = (de/dx), (x — x,) (12) 


where x,, it will be recalled, is the total chordwise extent 
of the roughness, and the subscript x, to de/dx denotes 
the value at the limit « = 0. Further, as remarked in 
Section (1), the linear part of the curve may be calcu- 
lated satisfactorily on the basis of certain assumptions 
relating to the nature of the particle impact. Also, the 
field of flow surrounding the airfoil is taken to be irrota- 
tional (a valid approximation at the small angles of in- 
cidence concerned), in which case (de/dx),, reduces to a 
function of the incidence and of certain terms associated 
with the geometry of the airfoil and the particle path 
at the point x = x,. Thus, ¢ can be determined, once 
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x, has been calculated, on the lines indicated in Section 
(1).* Alternatively, it will be clear, from the above re- 
marks, that model measurements of € will hold under 
full-scale conditions provided x, remains constant, 
though the reference in Section (3) to the somewhat un- 
representative nature of the relevant wind-tunnel data 
should not be overlooked.t| Nevertheless, since x, de- 
pends primarily on the ratio of the minimum rupture 
velocity and the free-stream velocity, it is not often 
that the range which can be covered in a wind tunnel 
coincides with that appropriate to the full-scale. This 
is the case regarding the present measurements in rela- 
tion to examples I, IT, III. 


Accordingly, since theoretical estimates of x, and « 
for the conditions I, II, III would have to be presented 
here without proof or experimental support, we shall 
resort to an empirical procedure which, though of a 
somewhat approximate nature, will be quite adequate 
for the immediate demonstration. Fig. 22 is a plot of 
[de/d(x/c)],, with respect to c/x, taken from the curves 
of Figs. 4, 5, 6, and 8. With the possible exception of 
the flat plate observation, the relation may be treated 
as linear to the order of accuracy which we require. 
Hence, if x, is known, e can be estimated from Eq. (12) 
and Fig. 22. For data in the former connection, we 
refer to reference 5, where the total extent of the rough- 
ness on both surfaces of the wings of several aircraft is 
recorded from observations in flight. Thus, in three in- 
stances for which the local thickness-chord ratio, t/c, 
of the wing varied between 0.090 and 0.095, the value of 
x, ranged from 0.30 to 0.39 of the chord,{ and for the 
thickest section investigated (¢/c = 0.175), the length 
x, was found to be 0.25c. This decrease of x, with in- 
crease of airfoil thickness accords, it will be noted, with 
the wind-tunnel records in Fig. 8 of the present paper. 
Strictly, for a given airfoil section, there will, in general, 
be three values of x, corresponding to the conditions I, 
II, III, and hence, three curves of «. However, to the 
order of accuracy which is sufficient in the present dis- 
cussion, the relevant differences can be ignored, and the 
above figures for x, may be taken to apply with equal 
validity in each case. Accordingly, the curves of €,;i: 
for the cases I, II, III were determined in relation to the 
NACA 66-009 section from the distribution of U’ given 
in Fig. 16 at an incidence of 1°, since this may be re- 
garded as about the attitude for cruising speed, and cor- 
respondingly, the envelope of € was estimated for the 
largest of the above values of x, appropriate to a section 
of this thickness, namely, when x, = 0.39c. 


* Such calculations depend on the value of « [see Section (4)], 
but the precise value, within the normal limits, does not affect 
the result appreciably. 

¢ Calculations show, however, that the error is small, the 
laboratory observations of x, being somewhat in excess of the 
corresponding free air values. 


t Strictly speaking, the measurements were made along the 
surface, but the difference between this distance and the corre- 
sponding projection along the chord is so small that it can be 
ignored. 
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The results presented in Fig. 23 show that the extent 
of the significant roughness, x,, amounts to about 0.06c 
in case I, 0.07¢ in case II, and 0.10c in case III. These 
figures agree well with Johnson’s predictions in refer- 
ence 5. Hence, even at the extreme of full-scale condi- 
tions, the amount of surface which has to be protected 
against roughness is quite small. 

The curve of € corresponding to the NACA 66-018 
section (x, = (0.25c) has also been plotted in Fig. 23, but 
the associated curves of €,,;,; have not been determined, 
since the small changes in LU’ due to the difference of 
section will yield results which deviate only slightly 
from those obtained for the NACA 66-009 section. In 
this case, it will be seen that the significant roughness is 
confined at most to some 7 per cent of the chord from 
the leading edge, and is rather less than the extreme 
value for the NACA 66-009 section due to the appre- 
ciable reduction of x, with increase of airfoil thickness. 

Finally, a curve of « has been plotted in Fig. 23 for 
the flat plate. Since no full-scale data in regard to x, 
are available for the case when ¢/c — 0), resource to a 
theoretical estimate of « has been necessary. This was 
based on the largest value of x,—namely, that relating 
to III in Table 2—for which the roughness was esti- 
mated to have spread over virtually the whole of the 
lower surface. Thus, the flat plate curve shown in 
Fig. 23 represents an upper bound, and indicates con- 
clusively that the very thin airfoil is almost completely 
free of significant roughness up to the highest Reyn- 
olds Numbers encountered in flight. 

There is also the interesting result, arising from these 
calculations, that the very thin and relatively thick sec- 
tions are less seriously affected by natural roughening 
than is the section of medium thickness. At the same 
time, the significant roughness is never more than a 
moderate percentage of the total roughness. 


(9) CoNCLUSIONS 


From the experiments described in Sections (4) and 
(6), and from the analyses of Sections (7) and (8), it is 
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Fic. 22. Empirical evaluation of « for natural roughness over 
that part of the envelope where the distribution is approximately 
linear. 
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Fic. 23. Examples of the degree of significant roughness at 
extreme values of full-scale Reynolds Number. 


TABLE 2 
Reynolds Approx. Wing 
Altitude, Number, Speed Lo, Chord, 
Case ft. R. X 10°* m.p.h. c ft. 
I 25,000 10 250 8 
Il 30,000 25 400 15 
III 30,000 50 600 20 


evident that, although the actual distribution of ex- 
crescence height in the case of roughness from in- 
sects is extremely irregular, the effect of such disturb- 
ances on the boundary layer may be consistently ex- 
pressed in terms of ¢ as defined by the envelope to the 
curve—that is, only the local extremes of ¢ are signifi- 
cant. Further, it appears that premature transition 
takes place in essentially the same way as it does for ar- 
tificial roughness of a three-dimensional character—.e., 
transition occurs at a roughness element when the local 
Reynolds Number X, attains a critical value which, for 
all practical purposes, depends only on the chordwise 
position of the element provided the chord Reynolds 
Number is not too small (in excess of 10’ approxi- 
mately), and the velocity distribution along the chord 
corresponds to that of a low-drag section at a small an- 
gle of incidence. 

However, in contrast to artificial excrescences such as 
rivets, lap joints, etc., for which « is likely to remain 
fairly constant over the surface, the distribution of 
natural roughness is uniquely determined, under given 
flight conditions, by the dynamic and aerodynamic re- 
lations associated with the airfoil and insects trapped on 
the surface. The resulting pattern takes a character- 
istic form—namely, one exhibiting a high concentration 
of roughness in a confined region near the leading edge, 
followed downstream by a relatively large area in 
which the deposits are much less densely packed. In 
terms of excrescence height ¢«, these two regions are 
sharply defined. A maximum is found to occur close to 
the forward stagnation point, but in a distance of | or 2 
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per cent of the chord further downstream ¢€ decreases by 
a matter of 60 to 80 per cent. From this position to the 
limit where ¢ is zero, the envelope of the high-spots is 
roughly linear. Thus, the possibility of a disturbance 
of the boundary layer is primarily determined by the 
leading-edge roughness, for which e may always be suf- 
ficiently large to cause transition. However, if the ef- 
fect of this initial region is suppressed by artificial 
means, as discussed in Section (6), then, even at the 
highest Reynolds Numbers occurring in flight [see 
Section (8)], most of the remaining roughness (of the 
order of 75 per cent of the total for an airfoil of typical 
thickness) is insignificant—i.e., the surface, though 
roughened, is aerodynamically smooth. This conclusion 
confirms the conjecture made in reference 5. 

It has also been found that the resultant velocity pro- 
file in the boundary layer immediately downstream of a 
roughened area depends appreciably on the intensity of 
roughness (proportion of unit area of surface covered by 
deposits). Thus, if the critical deposits are widely 
spaced, the diffusion of the individual wakes is not suf- 
ficient to produce a uniformly distributed boundary- 
layer profile characteristic of fully developed turbu- 
lence, but leads to one intermediate between this state 
and that corresponding to purely laminar flow. For 
example, at an incidence appropriate to cruise or climb 
conditions, the first signs of general departure from the 
laminar profile have been observed to take place when 
the roughness has spread to about 2 per cent of the chord 
from the leading edge, and when the maximum inten- 
sity is approximately 5 per cent, but by the time the 
leading-edge intensity has risen to some 30 or 40 per cent 
and the chordwise extent of the roughness is rather more 
than half the maximum attainable under the relevant 
flight conditions, the fully developed turbulent bound- 
ary layer may be expected to exist everywhere in the re- 
gion downstream of the roughness. This effect of 
roughness intensity on the ultimate state of the bound- 
ary layer, once mixing of the wakes from individual ex- 
crescences has taken place, is not a factor, of course, in 
the promotion of transition which, as we have observed 
above, depends only on the local roughness Reynolds 
Number X,. 

With the possible exception of a small interval in the 
region of zero lift, the total chordwise extent of the 
roughness (x,) at a constant wind speed is an approxi- 
mately linear function of incidence over the range ap- 
propriate to normal flight conditions. Theory leads to 
the same result, and shows that it is a fair approxima- 
tion even for the speed-incidence relation peculiar to a 
given aircraft within the limits of maximum speed and 
take-off or landing. Thus, since the rate of change of x, 
with incidence and speed is positive for the lower sur- 
face, the maximum spread during flight, under condi- 
tions favorable to the promotion of natural roughening, 
is determined approximately by the maximum incidence 
and minimum speed at which the aircraft is flown. On 
the other hand, the opposite holds for the upper sur- 
face—i.e., x, decreases with increase of incidence and 
reduction of speed—so that the extent of the roughness 


then assumes a maximum in the region of the maximum 
speed and least incidence reached during a flight. 
Finally, it appears, though the evidence is at present 
very limited, that the wide differences in the physical 
and mechanical properties of insects encountered in 
flight have only a secondary effect on the distribution 
and magnitude of ¢ over that part of the surface where 
the roughness is aerodynamically insignificant. This is 
not unexpected, since deposits in the region concerned 
are composed of drops of body fluid and small visceral 
particles. Near the leading edge, however, where 
whole insects are found to adhere, this conclusion no 
longer holds, but, as we have remarked earlier, the size 
of the excrescence in this area is always so large, even 
for the very small and fragile insect, that transition to 
turbulence will almost certainly occur close to the lead- 
ing edge if artificial means are not used to circumvent 
this premature collapse of the laminar boundary layer. 
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Chemically Frozen Boundary Layers With 
Catalytic Surface Reaction 


DANIEL E. ROSNER* 


Princeton Unwersity 


SUMMARY 


When a fluid boundary layer containing a single reactant de- 
velops along a solid surface which acts chemically as a sink dis- 
tribution for this species, its steady-state concentration along 
the surface, and, hence, the reaction rate distribution depends 
both upon the true surface chemical kinetics and a combination 
of convective and diffusive mass transport through the fluid. 
Considered here is a simple approximate method for studying this 
streamwise distribution of reaction rate in terms of the shear and 
temperature distributions along the surface. The method is 
based upon a treatment of the concentration integral equation 
which initially parallels Liepmann’s recent treatment! of the 
energy integral equation but which ultimately leads to a differen- 
tial equation for the concentration distribution of reactant along 
the surface. This approach to ‘‘mass transfer’’ problems, which 
is intended to bring out clearly the role of the governing catalytic 
parameter as well as several interesting properties of solutions, 
is then compared with alternate techniques, in particular with 
that of Chambré and Acrivos.2~4 The reactant concentration 
is considered to be so small that constant diffusion coefficient, D, 
viscosity coefficient, 4, Prandtl Number for diffusion (Schmidt 
Number), Prp = u/(pD), can be assumed, and the fluid is taken 
to be incompressible, with catalyzed reaction taking place only 
at the wall. Thus, the hydrodynamic field is unaltered by the 
surface reaction, and only the effect of small surface temperature 
differences on the kinetics of the wall reaction—i.e., the effective 
“sink”? strength—are taken into account. Extension to the 
case of several reactants is discussed, and an interpretation of the 
governing catalytic parameter as a ratio of characteristics times® 
is given. 


(1) INTRODUCTION 


i 1950, LIGHTHILL® OBTAINED an expression for the 
normal temperature gradient distribution along a 
surface in terms of the surface shear and temperature 
distributions which is asymptotically exact in the limit- 
ing case of infinite Prandtl Number. The method is 
based upon an operational treatment of the boundary- 
layer heat convection equation in the von Mises form, 
into which a local approximation for the velocity field 
is introduced. If, along with the shear stress, the sur- 
face temperature is considered to be known every- 
where, the heat-transfer distribution is thereby deter- 
mined. The equation for convection and diffusion of 
mass in a boundary layer, being the same as that for 
heat convection, therefore yields to the same attack— 
i.e., if the hydrodynamic boundary condition at the 
wall is unaltered by the mass transfer, then Lighthill’s 
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result may be applied to the problem of surface reac- 
tions, as was pointed out by Chambré.* The difference 
in the classical statement of the two problems is that, 
while one usually finds the heat-transfer distribution 
corresponding to a prescribed surface temperature dis- 
tribution, for the case of surface reactions the concen- 
tration distribution of reactant along the surface is the 
essential unknown, and the reaction rate at the surface 
is itself presumed to be independently expressible, 
according to ‘“‘heterogeneous’’ chemical kinetics, in 
terms of the reactant concentration at the surface. 
There results an integra] equation for the reactant con- 
centration, the solution of which is discussed in refer- 
ences 2 and 3. This is analogous, therefore, to the 
calculation of the steady-state temperature distribution 
of freely radiating solids.® 

Using integral methods, Liepmann'! and Chambré‘ 
have recently given simple derivations of ‘Lighthill- 
type” formulas which are of pedagogical interest in 
that the physical features of the problem are not ob- 
scured by details of the operational method. Chambré 
obtains a result [Eq. (1S8b)] using the von Karman- 
Pohlhausen integral method (with cubic profiles) which 
is stated to agree with Lighthill’s formula [Eq. (18a) } to 
within 1 '/. per cent. In reference 3, use is made of 
the Lighthill relation itself in order to obtain an integral 
equation for the reactant concentration and, hence, 
reaction rate distribution along the surface. However, 
apart from basic differences in method, the assumption 
(3) underlying the Lighthill result should be distin- 
guished from the assumptions (2) and (9) which underlie 
Liepmann’s recent treatment.' In this paper, Liep- 
mann’s approach is adapted to the problem of surface 
reactions since, with the aid of the quasi-similarity as- 
sumption following Eq. (8), the resulting mass transfer 
equation may be treated as a first-order differential 
equation to better illustrate the behavior of solutions. 
The emphasis here will be upon the display and inter- 
pretation, in a simple manner, of the role of the relevant 
transport and catalytic parameters. This is to provide 
at least a qualitative basis for the understanding of more 
complex problems of this general class. 


(2) APPROXIMATE SOLUTIONS 


The starting point is the concentration integral 
equation for the two-dimensional diffusion boundary 
layer, which may be written 


4 
(d/dx) pu(c, ody! (1) 
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Here R,, is the rate of surface consumption of the react- 
ant which has the local concentration (weight fraction) 
c within the layer and c, in the free stream. If the 
Prandtl Number, ?7p, for diffusion is not very small, 
then the reaction rate, R,,, at the wall should be prin- 
cipally determined by the velocity profile in the im- 
mediate neighborhood of the wall, and not by its asymp- 
totic behavior; as was first suggested by Fage and 
Falkner’ (1931), later used by Lighthill (1950) for the 
case of heat transfer, and recently adapted by Chambré 
and Acrivos (1956-1957) to the problem discussed here. 
With this in mind, for laminar flow far from separation 
we express the local velocity u in terms of the local con- 
centration c and the shear stress 7,,, using c as a Crocco 
variable: 


u(c) = (€ — Ce) (Prp) (Ru)} (2) 


That this is valid very near the wall is a consequence of 
the local relations 


u(y) (Tw/u)y (3) 
and 
c(y) =~ Cy + (Oc OV) = Cw + {Ro (Dp) y (4) 


although it is significant that the elimination of more 
general functions of the ordinate, y, can also lead to 
Eq. (2). Therefore, the relation (2), and not (3) and 
(4) individually, should be taken to comprise an under- 
lying assumption of the method. 

If we introduce the local flux, R, defined by 


R = Dp(0dc/dy) (5) 


(the magnitude of which is R, when evaluated at the 
wall), then the integration indicated in Eq. (1) can be 
performed over c(y) to give 


R,, = (Dp?/Prp) (d/dx) X 
Ref — Cy) (Ce — ©) (Rr, Rack (6) 


With the further introduction of a “‘similarity’’ vari- 
able n(c), and a streamwise variable ¢(x), defined, re- 
spectively, by the relations: 


n(c) = (¢ — Cw)/(Co — = (7) 


(where Z is a characteristic body dimension) Eq. (6) 
becomes: 


Ry = |Dp?/(2LPrp)} (1/)(d/de) X 
1 
~ — n) (Ry Rant (8) 


If it is tentatively assumed that R/R,, is a function of 
n alone, then the integral in Eq. (8) becomes a pure 
number, which will be denoted by ap. The conse- 
quences of this ‘‘quasi-similarity’”’ assumption and the 
evaluation of a» will be examined later, in the discussion 
of a special case. There results the basic equation: 


Ry = | aDp?/(2LPrp)} (1/g)(d/de) X 
{ Rr} (9) 


At this point the kinetics of the surface reaction is 
introduced with the postulate that the reaction rate R,, 
at the wall is independently expressible as a function 
of the local partial density pc, of reactant. Eq. (9) 
may then either be regarded as a first-order nonlinear 
differential equation or treated as a nonlinear integral 
equation for the reactant concentration c,(¢) along the 
surface. The equation is general, with respect to stream- 
wise changes in the catalytic condition of the surface 
and free-stream concentration, and the limitations 
upon its utility are the extent to which the streamwise 
gradients, appearing in the solutions, are in accord with 
the underlying assumptions of boundary-layer theory, 
and the extent to which both Eq. (2) and the kind of 
“similarity” postulated in obtaining Eq. (9) are rea- 
sonable approximations. 

The integral equation for the normalized reactant 
concentration, = ¢»/c., is here expressed in terms 
of the “‘Nusselt Number, Vip, for mass transfer’’ for 
the case c, = constant: 


Nup = (0.332a/3)'* (Prp)'* (Re,)'’? X 
Vir = f = 


-1/3 


(10) 


where S(¢) is a dimensionless ‘‘shear function’’ defined 
by: 


Tu(f) = 0.332 (11a) 


Nup = R, |Dple. — Cw), /x} (11b) 


and Re, is the Reynolds Number based on the length, 
L, and approach velocity, U’. 

There is a comparable result for the case of pure heat- 
transfer with constant density-viscosity product, pu. 
If h is the normalized surface thermal enthalpy, h,,/h., 
and ay is the heat flux analog of the ay defined above, 
then Liepmann’s result may be written in the dimension- 
less form: 


Nu, = (0.332 (Pr)? (Re,)'? X 
—1 
— Des) eV (1 = ac’ 
0 
(12) 


where Nz, is the ordinary local Nusselt Number for 
heat transfer, Pr, being the Prandtl Number for heat 
conduction. If the a’s are identified, then it is ob- 
served that the Nusselt Number for mass transfer may 
be obtained from the heat-transfer Nusselt Number by 
formally making the replacements: Pr, > Prp, h > 


This is, in fact, the basis for the simplest possible 
approach to this class of problems, as is briefly discussed 
later. One can imagine cases where both thermal 
energy and mass are transferred, and for which Eqs. 
(10) and (12) are valid simultaneously. Moreover, if 
the reaction at the surface is exothermic, then the net 
energy transfer will be the sum of contributions from 
local molecular conduction, and the diffusion of reac- 
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Fic. 1. Integral curves of the catalytic flat plate equation. [Eq. 
(16) with S = 


tive species to the wall. If the surface temperature 
distribution were prescribed, then the corresponding 
thermal enthalpy distribution could be readily esti- 
mated, and the conductive heat-transfer distribution 
would follow from Eq. (12) directly. However, Eq. 
(10) would remain an integral equation in content, 
since the steady-state distribution of reactant along 
the surface itself depends upon the kinetic law Ri(Cx)- 
Thus, in order to determine the contribution of surface 
reaction to the local energy transport, one must first 
find the concentration distribution c,(¢), and then 
introduce this result into the kinetic law. 

In order to display the nature of solutions and the 
role of surface chemistry, Eq. (9) will be treated in- 
stead as an ordinary differential equation, and the fol- 
lowing explicit reaction rate law will be assumed: 


Ry = ky { pC} ” (13) 


Here k, represents a wall temperature dependent 
“catalytic rate constant’ and n is referred to as the 
“order’’ of the chemical process. Specification of k,, 
and » thus completely describes the chemistry of the 
problem. For the purposes of illustration, k,. will be 
considered constant. If one then defines a ‘‘stretched”’ 
streamwise variable z by: 


z = (14) 
where @,, is the composite ‘‘catalytic parameter” 


C,, = { (D/L) (Re,)'?}-! (15) 
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it is found from Eq. (9) that the normalized reactant 
concentration ¢ = c,/c, satisfies the first-order, ordinary 
differential equation: 


dé/dz = {[—@(1 — @)]/[2n + &(3 — + 
(22/S) [e"/(1 — — (1/S) (dS/dz)! (16) 


where: S(z) = S[¢(z)]. 

The “catalytic parameter” appearing implicitly 
in Eq. (16), is fundamental to this class of problems. 
It is observed from this equation that flows for which 
S is a constant have the similarity property so that the 
normalized reactant concentration € depends upon 
€, and ¢ only through their product @,¢. This point 
will be returned to in the discussion of the catalytic 
parameter. 

The diffusion of the dilute reactant and its catalyzed 
reaction at the surface will not be considered to alter 
momentum transfer in the layer appreciably. There- 
fore, S(¢) will be independently determinable from 
approximate analytical procedures (such as that due 
to von Karman and Pohlhausen), from exact solutions 
of the boundary-layer equations, or from hydrodynamic 
experiments for the related nonreactive system. For 
example. if one considers the laminar boundary layer 
which develops from the leading edge of a semi-infinite 
catalytic flat plate, then S is simply unity, which, of 
course, prompted the definition (11a). 


(3) QUALITATIVE BEHAVIOR OF SOLUTIONS 


To illustrate the solutions of Eq. (16) in a particular 
case a further specialization to the case of the semi- 
infinite catalytic flat plate is made. Then the normal- 
ized concentration distribution ¢(z) along the surface 
satisfies Eq. (16) with S set equal to unity. Concern- 
ing the initial condition on ¢, we specify a priori that 
for bodies with sharp leading edges ¢(0) = 1. This is 
equivalent to the statement that, regardless of reaction 
order, the concentration boundary layer starts with 
zero thickness. It is significant for what follows to 
note that the ‘‘initial’’ point (1, 0) of the és plane is 
singular. In Fig. 1 is shown the qualitative behavior 
of the integral curves for arbitrary reaction orders 
n greater than 3/2. For positive reaction orders 
smaller than this value there is no saddle point in the 
first quadrant of the é-z plane; the qualitative behavior 
of integral curves in the domain of physical interest: 
1 > € > 0; z > O, however, is unchanged. As de- 
picted, there exists a one-parameter family of integral 
curves passing through the initial point, and these are 
locally nonanalytic. However, a distinguished curve 
emerges from this point which 7s analytic, its slope 
there being —(1/2)'* for all reaction orders. If, 
then, there is only one solution to the problem posed, 
we ask which of these eligible curves in fact is the solu- 
tion. To decide this, all candidates locally incon- 
sistent with one or more of the assumptions underlying 
the method must be discarded. From Eq. (2), the 
only integral curve for which the velocity ~ near the 
plate has a finite limit as z ~ 0 is that for which | — 
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Fic. 2. Normalized concentration distribution for uniform 


catalytic plate with first-order surface kinetics. A comparison 
of three methods. 
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€ ~ s. On this basis, the analytic curve may be 
elected as the solution to the uniform catalytic plate 
problem. 

It is now interesting to consider more general prob- 
lems for which these “extraneous” integral curves have 
physical meaning. An example would be the ‘“‘piece- 
wise-uniform” catalytic plate—i.e., a plate which, at 
distinguished values ¢,*, ¢2* . . . of ¢ is divided into 
distinct regions of different catalytic activity. Then, 
at each critical value z* of z one would jump on to the 
integral curve which passes through the same value of 
€ but at the new value of z. Thus, integral curves lying 
to the right and left of the analytic solution become 
physically meaningful when the plate increases or de- 
creases its catalytic activity at each z*. Furthermore, 
if one considers as the limiting case of a ‘‘piecewise- 
uniform’’ catalytic plate a surface (say, nonisothermal) 
for which the catalytic activity varies continuously, 
then the solution for this concentration distribution 
would correspond to a continuous sequence of jumps 
on to neighboring integral curves. 


(4) A QUANTITATIVE COMPARISON OF METHODS 


In order to compare a result of this method with a 
numerical result of alternative methods, it is necessary 
to select a reasonable value of ao [see Eq. (8) ] which, in 
turn, depends upon the choice of the function R/R, = 
g(n). An extreme choice: g(n) = 1 — 7 leads to the 
value a) = 1/2; a more realistic choice g(n) = 
/(1 — 7?) leads to ay = 0.215. In any event, since 
a is raised to the 1/3 power, the choice is not critical 
regarding the relation of € to z. The analytic solution 
curve of Eq. (16) for the first-order case, obtained by 
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forward numerical integration, is compared in Fig. 2 
to the analytic solution of Chambré and Acrivos, and 
an “engineering”’ solution, using the value aj) = 0.215. 
The curve marked ‘‘Eq. (17),” with ¢(0) = 1, is thus the 
analytic solution of the equation 


dé/dz = {—e(1 — + 2)} x 
{(1/z) + — (17) 


As already indicated, Chambré and Acrivos have 
employed the mass-transfer analog of the Lighthill 
heat-transfer equation. In the notation adopted 
here, this equation may be written: 


Nup = 0.296(Pro)' (Re,)"? — X 
— -—1/3 


With the initial condition ¢(0) = 1, this integral equa- 
tion yields the desired reactant concentration distribu- 
tion if the surface shear} and reaction rate law are 
specified. It is stated to be sufficiently accurate 
for values of Prp greater than or equal to 0.7 and, like 
Eq. (9), holds for nonseparating flows. In reference 
2, these authors have outlined a numerical technique 
for constructing solutions to this integral equation for 
arbitrary shear distributions and chemical kinetic rate 
laws. However, for the first-order, isothermal plate 
case a convergent power series for the analytic solu- 
tion ¢(z) is given,* which is shown in Fig. 2. It is 
interesting to note that Chambré‘ has obtained a mass- 
transfer equation very similar to Eq. (18a) by using the 
von Karman-Pohlhausen integral method with cubic 
profiles. In the dimensionless notation adopted here, 
Chambré’s integral equation may be written: 

Nup = 0.292(Prp)"* (Rex)"” [g2/(1 — @] X 

t t } -1/3 


This should be compared with both Eqs. (12) and 
(18a). Only Eq. (18a) has been employed* to obtain 
the concentration distribution along the surface of a 
catalytic flat plate. 

Finally, these methods are compared to an ultra- 
simplified approach to this class of problems, having 
its origin in the so-called ‘“‘stagnant film’ hypothesis. 
Here a stagnant film is imagined, representing the entire 
resistance to transfer of energy or mass to the surface. 
Accordingly, if one knows the effective film thickness 
“6” from a heat-transfer solution or experiment, then 
the mass-transfer rate under ‘‘comparable’’ conditions 
can be obtained by identifying film thicknesses, pro- 
vided the roles of the Prandtl Numbers Pr, and Prp 
are interchanged when they are not equal. Identifi- 


cation of the film thicknesses gives 
MT. — Tw)/g, = = Dolce — Cy)/m" (19) 


Actually, following Tifford,! Chambré and Acrivos suggest 
the use of an effective shear 7, related to the actual shear tw by 
the local pressure gradient and momentum thickness of the 
boundary layer. 
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CHEMICALLY FROZEN 


where m” represents the local rate of mass transport 
per unit time and area, and g, represents the local 
heat-transfer rate; \ being the thermal conductivity of 
the fluid. In terms of the corresponding local Nusselt 
Number Nuy, Eq. (19) then gives, for the mass transfer 
rate: 


= Nu,(D/L)pc.(1 — (20) 


The laminar, flat-plate Heat transfer problem with con- 
stant wall temperature has been studied analytically 
by Pohlhausen® (1921), who obtained the “self-similar” 
temperature profiles corresponding to the Blasius ve- 
locity field. His numerical results could well be ex- 
pressed by the relation 


Nu = 0.332 (Pr,)'? (Re,)'*¢ (21) 


Since Eq. (20) is presumed to hold when steady-state 
surface reaction occurs, then m” must be identified 
with the rate R,, of reactant consumption at the wall. 
If we now disregard the nonsimilarity correction at- 
tributable to the surface boundary condition, then 
the Nusselt Number Nu, given by Eq. (21), may be 
introduced into Eq. (20), upon replacing Pr, by Prp. 
The normalized reactant concentration ¢ may then be 
obtained from the algebraic equation 


e,¢é" — 0.332}1 — = 0 (22) 


The result for m = 1 is plotted in Fig. 1. It may be 
remarked that, provided c, Cy 18s chosen as the 
“driving force’ for mass transport, an identical result 
is obtained if at the outset one alternatively draws 
upon the analogy between momentum and mass trans- 
fer in fully developed turbulent flow, and the Chilton- 
Colburn modification of this analogy to account for 
molecular transport processes. This modified analogy 
may be stated in the symmetrical form:' 


t,/(v?3U) = m"/[D? (ce. — Cy) ] (23) 


which provides a simple relation between m” and the 
wall shear stress. If Eq. (23) is assumed to hold 
locally, then introduction of the Blasius shear distribu- 
tion also gives Eq. (22). For more general shear dis- 
tributions, S(¢), the analogy provides us with the rela- 
tion: 


e,¢é" — 0.332S(¢) {1 — = 0 (24) 


solutions of which are easily expressed for n = 1, 2. 
Each of the above solutions has the property that, 
at the leading edge, where the shear stress is unbounded, 
the reaction rate is ‘‘chemically controlled,” R,, there 
being Rut pce}, regardless of diffusion or convection 
parameters. Far downstream, however, where both 
the shear stress and reactant concentration at the wall 
approach zero, the reaction rate may be said to become 
“diffusion controlled’”’ in the sense that the catalytic 
activity of the surface no longer enters the result— 
i.e., the surface has the kinetic capacity for greater 
conversion but must content itself with reacting only 
what can reach it from the free stream. With regard 
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to the ‘‘apparent’’ kinetics, the reaction rate near the 
leading edge therefore responds to changes in free- 
stream reactant concentration in a way which reflects 
the true m order surface kinetics. Further downstream 
the diffusion process masks the true reaction order. 
Ultimately, the reaction rate responds to changes in 
free-stream reactant concentration as if a first-order 
reaction were taking place, regardless of the true reaction 
order. Upon comparison of these three approaches, it 
is seen from Fig. 2 that the numerical differences are not 
large for this particular example. In other cases, the 
choice of method is more critical,’ and modified integral 
techniques of the type discussed here offer alternate 
methods of both intermediate complexity and accuracy, 
capable of generalization to include turbulent flows. t 

Extension of these treatments to the case of more 
than one reactant is made by considering each reactant 
concentration to be governed by its own Eq. (9) with 
kinetic rate, diffusion coefficient etc., corresponding to 
the species in question. For instance. if species A and 
B present in the free stream react at the surface accord- 
ing toA + B> Cwith Ra w = Rew = e*kw\caca} w 
and each free-stream concentration is specified, then 
the concentration distributions c4 w(t), cg and, 
hence, the rate Ry(¢) may be determined. In the 
special case when a»D* is the same for both components, 
and the free-stream concentrations are equal as well, 
the normalized wall concentration ¢(¢) of each species 
satisfies Eq. (16) with nm = 2. 


(5) THe CaTALyTIC PARAMETER AS A RATIO OF 
CHARACTERISTIC TIMES 


The catalytic parameter C,, is seen to enter this class 
of problems in a natural way. As discussed by Chambré 
and Acrivos,* the magnitude of this type of parameter 
determines the streamwise extent of the effective 
“chemically controlled” region of a developing flow. 
This is displayed in Eq. (17) by the fact that the nor- 
malized reactant concentration € = c,,/c, is a function 
of the product @,¢ alone. Thus, if the catalytic 
parameter @,, is large, a given € < | is attained at a 
fractional length x/Z which is correspondingly less 
than the value for a small @,,. 

We indicate here that this ‘catalytic parameter’ 
may be looked upon as a ratio of characteristic times— 
i.e., the ratio of a time raire, for the reactant to diffuse 
across the concentration boundary layer of thickness 
6 at L, to a time r"chem. representative of the rate at 
which the wa// can consume a prescribed amount of this 


t Added in proof: Ambrok’s approximate method (Sov. Phys., 
Tech. Phys., Vol. 2, No. 7, A.I.P. Translation, September, 1957, 
pp. 1979-1987) has been found to be particularly convenient for 
this class of problems. Its application to the catalytic plate 
leads to closed form results for all reaction orders. (Rosner, 
D. E., AeroChem Technical Memorandum No. 12, December, 
1958.) In terms of the variable z'(+/2/0.332)@,¢ the correspond- 
ing solution for first-order kinetics is 


(2’)?/2 = [(1 — &)/e*] — In (1/e?) 


If plotted on Fig. 2 this solution would be almost coincident with 
the power series solution given in reference 3. 
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reactant.’ In this respect, the catalytic parameter may 
be regarded as a kind of ‘‘wall Damkoéhler Number.” 
For example, consider the plate problem, and for 
7" chem. Choose a time proportional to the time that it 
would take a unit area of the surface to consume the 
reactant lying immediately above it in the boundary 
layer. Then: 


chem. ~ 5k, '(pc,)~ "~? (25) 


Since the time to diffuse across the layer goes as 6°/D, 
and, as indicated by Eq. (21), 6/L ~(Prp)~ 
we indeed find that €,, is a measure of the ratio of these 
times. The catalytic parameter enters more com- 
plicated problems of this general type but somewhat 
in disguise. For instance, the group 


RuPu( Pro)? { peme(du,/dx) ro} (26) 


which enters recent hypersonic compressible stagna- 
tion-point solutions'!: |? is essentially the catalytic 
parameter @,(for first-order surface reactions) based 
on the nose radius R, and the fluid velocity U’; behind 
the bow shock. This can easily be verified in the limit 
p = constant, with the inviscid velocity gradient at the 
nose going as U,/R,. These remarks bear upon the 
question of scaling or simulation. To the usual simili- 
tude requirements for aerothermochemical systems'* 
one must, in general, add a dimensionless “‘criterion”’ 
like the catalytic parameter C,,—.e., a ‘‘wall Damkohler 
number’ to account for the role of ‘‘heterogeneous”’ 
chemical reactions in the physics of such problems. 
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Note on a Criterion for Severity of 


Roll-Induced Instability 


ORDWAY B. GATES* ann KARLIS MINKA** 
The Marlin Company 


SUMMARY 


A method is presented which may be used to evaluate the 
severity of the roll-induced instability inherent in many current 
aircraft designs. The analysis is based on the assumption of a 
steady rolling velocity, and for this case, expressions are pre- 
sented from which the divergence rates can be estimated. Also, 
from knowledge of the divergence rates, the magnitude of the 
divergence may be estimated as a function of the angle through 
which the aircraft is rolled. 


SYMBOLS 


m = mass, slugs 
I,,l,, 72 = principal moments of inertia, slug-ft.* 
’,4,r” = angular velocities about x, y, and s axes, rad./sec. 
a = angle of attack of principal x axis, rad. 
B = angle of sideslip, rad. 
L, M, N = rolling, pitching, and yawing moments, ft.-Ib. 
1.2 = forces along y and z axes, lb. 
Lg = 
= OL/op 
= OL/dr 
Lia = 
Me = 0M/da 
M, = 0M/0dq 
N, = ON/odr 
Neg = 
Vp = OY/0p 
= 02/d0a 
> = bank angle, rad. 
d = characteristic root 
s = Laplace transform operator 
t = time, sec. 
ba = total aileron deflection, rad. 
Subscripts 
0 = initial value 
ss = steady state 
©, C2 = denote critical values of roll rate 


A dot over a quantity indicates differentiation with respect to 
time. 


INTRODUCTION 


[ HAS BEEN DEMONSTRATED theoretically and verified 
in flight tests that, for certain combinations of roll- 
ing velocity, inertia characteristics, pitch-yaw damping, 
and pitch-yaw frequency relationships, aircraft can 
exhibit dangerous divergent (unstable) tendencies in 
rolling maneuvers. Phillips! presents a method for 
predicting the range of roll rates for which a given air- 
craft is unstable, which is based on the assumption of 
a steady rolling velocity. Numerous papers have 
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attacked the problem from this point of view, whereas 
other authors* * have attempted to discuss it from the 
standpoint of steady-state response or autorotational 
states of motion. 

The purpose of this paper is to present a criterion for 
severity of roll-induced instability which is based on (1) 
rate of divergence and (2) the magnitude of the diver- 
gence as a function of the angle through which the 
aircraft is rolled. 


ANALYSIS AND DISCUSSION 


The following equations, from which the aircraft 
lateral-longitudinal responses to aileron inputs can be 
computed, result from the assumptions that the effects 
of gravity can be neglected and that the component of 
the linear velocity along the principal x axis is constant: 


q — pB + (Z,/mv)Aa = & 
p (ao + Aa) — r + (Vg/mv)8 = 8 
— pr + (M,/1,)Aa + (M,/T,)q = 4 
— + (N,/I)r + = 
+ (Lp/Iz)p + + = p 


These equations are written relative to principal 
axes, and the parameter a» is the inclination of the x 
axis to the flight path at ¢ = 0. From the first four 
of these equations, for the case where é = B = g = 
* = p = O, expressions can be derived for the steady- 
state responses 6, Aa, g, and r per unit ap in terms of the 
rolling velocity p. Typical curves are shown in Fig. 
1. Substitution of these values into the roll equation 
yields p,, in terms of 6,,, for given values of ap, and 
typical results are shown in Fig. 2. From Fig. | it can 
be seen that for p = 1.2 rad./sec. and p = 2.7 rad.. 
sec., the steady-state responses Aa, 6, r, and g tend 
toward infinity. These values of p, of course, have 
significance only for the case being considered, and in 
no sense should they be construed to represent any 
given aircraft. 

These values of p correspond to the values between 
which the Phillips steady rolling analysis predicts 
instability for the case being considered. 

If the steady-state values of Aa, 8, g, 7, p, and 6, are 
considered as possible dynamic trim conditions, a 
stability analysis yields the fact that those combin- 
ations which lie between p = 1.2 and p = 2.7 are un- 
stable trim combinations and hence do not actually 
exist, whereas the remaining combinations are stable 
trim points. 
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Steady-state relationship between roll rate p and]Aa. 
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CRITERION FOR SEVERITY OF ROLL-INDUCED 


Although the steady-state responses of Aa, 8B, gq, 
and r theoretically become very large as p tends toward 
1.2 rad./sec. (and infinite for 1.2 < p < 2.7 rad./sec.), 
the rate at which the responses approach these large 
values is certainly the more important factor. If it 
is assumed that the roll rate p is maintained at a con- 
stant value, the solutions of equations (1), as p > 1.2, 
take the form, in the transformed state, of 


xi(s) = aog(s)/s(s — A)h(s) 

where g(s) and A(s) are polynominals in s, and A, is the 
real characteristic root which becomes unstable for 
1.2< p< 2.7. The inverse transform of this function 
iS 
xi(t)/ao = + 

...= Ge" +... 
For P = Periticar, 4 = 0 and 

x;(s) = aog(s)/s*h(s) 

Ci + CG’t+... 


x,(t) ‘ao = 


If the terms shown are assumed to be the ones which 
contribute mainly to the system responses in the 
vicinity Of Persicar and in the unstable roll rate range, 
then for p < p,,, or for p., << p< Pa, 

%(t)/ao ~ Je = AC, 
and, for p = pz, 
C,’ 


%,(t)/ ao 


The coefficient Cs’ is determined from an evaluation of 
the residue at the pole of order 2 at s = 0 of the function 


x(s)/ao = g(s)/s*h(s) 
= g(0)/h(0) 


Consider the case first where p = p,,, and hence A; = 
0. The change in x; subsequent to ¢ = 0 is given by 


and is 


Ax,(t)/ao = [xi(t)/ao] — [x(0)/a0] = C,’t 
For 
(e 1) 


Pe, < p < Press Ax ;(t)/ ao = 


Rather than plot Ax; as a function of ¢, let 
t= o/p 
and, therefore, 


(Ax;/ao) (6/p) = C2'(o/p), = Pe, 


or 
(Ax,/ ao) (¢/p) 1), Pa < Pe 


These functions for Aa/ay and 8/ao are plotted in 
Fig. 3 versus the bank angle ¢ for various values of p. 
Similarly, the divergence rates a/ ap and B/a are plotted 
in Fig. 4 versus roll rate p. Also shown in Fig. 4 is a 
plot of A; versus p. It is apparent from these Figures 
that for the case considered the divergence of Aa is 
greater for the lower critical roll rate and the divergence 
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Fic. 5. Aircraft responses to aileron deflection, 360° roll. 


of 6 is much higher at the higher critical roll rate. 
This result is due to the fact that the directional static 
stability is higher than the longitudinal static stability 
for the aircraft considered here, and the roll rate reso- 
nates first with the longitudinal frequency. These 
trends would be reversed if the opposite were true. 

The main point to be made here is that the large 
steady-state responses associated with roll rates ap- 
proaching the unstable range and the theoretically 
infinite responses within the unstable region mean 
very little unless the degree of instability and rate of 
divergence is known. It is believed that the construc- 
tion of plots similar to Figs. 3 and 4 for a particular 
aircraft will afford a good indication of the severity 
of the roll coupling problem for that particular aircraft. 

One further point to be made concerns the curve of 
Pss versus 6,,, in Fig. 2. This curve indicates that 
as Pp > pz, ba, > ©. This result is very unrealistic, 
as shown by the following analysis. The value of 6,,, 
was obtained from the equation 


It is apparent that unless the rolling is allowed to con- 
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tinue until r and 8 approach their theoretically large 
steady values, 6,,, for all practical purposes will be 
primarily a function of p,, only, and will certainly not 
tend toward infinity. This is demonstrated by the 
analog computed response for 6,,,,, = 0.40 rad. for the 
aircraft and flight condition being considered. The 
responses in 8, Aa, and p and the input 6, are plotted 
in Fig. 5. The response in roll rate is seen to be essen- 
tially constant at p = 1.6 rad./sec. The responses 
in 8 and Aa are divergent for this roll rate, but the 
degree of instability and rate of divergence (particu- 
larly the 6 response) is not high enough to have a signi- 
ficant effect on the roll rate response even though the 
aircraft rolled through one complete revolution. 


CONCLUDING REMARKS 


A method has been derived which may be used to 
determine the severity of roll-rate-induced divergence 
inherent in present-day aircraft. It is pointed out that 
the intolerability of this type of instability is very 
strongly dependent upon the rates of divergence which 
are encountered. It is further shown that the diver- 
gence can occur initially in either the pitch or yaw mode, 
dependent upon the relationship between the natural 
frequencies of these modes. Expressions are presented 
which can be used to estimate the magnitude of the 
divergence of a particular mode as a function of bank 
angle. 
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Flight Mechanics of Low-Thrust Spacecraft 


FRANK M. PERKINS * 


Convair Astronautics, A Diviston of General Dynamics Corporation 


SUMMARY 


Spacecraft starting from circular orbits and using a low thrust, 
tangentially directed, will follow a spiraled flight path as illus- 
trated in Fig. 1. Whether or not the thrust is constant, if it re- 
mains less than one per cent of local vehicle weight the mean path 
of these trajectories will maintain the familiar circular velocity- 
altitude relationship for some time. Complete trajectory equa- 
tions for this mean path and approximate equations for the os- 
cillations about it are derived and presented here. 

For spacecraft using constant thrust acceleration, parameters 
are presented for generalizing all differential equations of motion 
to make them independent of the thrust acceleration and the 
gravity constant (choice of planet). In the case where this thrust 
is directed tangentially throughout flight, generalized plots are 
presented relating velocity, radius distance, and time. 

When the constant tangential thrust acceleration is less than 
10~? initial g, single curves independent of the value of the thrust 
acceleration result, relating the parametric velocity, radius dis- 
tance, and time. Plots of mean path velocity, distance from the 
center of the gravity source, and time may be constructed for any 
space craft of this type regardless of what planet, what initial 
altitude, and what thrust acceleration is employed, by simply 
multiplying the scales of the parametric plots by suitable con- 
stants. 

Simple algebraic expressions for speed, altitude, and time to 
achieve any specific point on the trajectory are available from the 
generalized plots. Examples of these are presented for the para- 
bolic escape velocity point. 


SYMBOLS 


a = dimensionless thrust acceleration = f/(u/ro?) 
f = thrust acceleration = F/ 
F = thrust 
g = gravity acceleration (at any altitude) 
h = altitude above surface of planet 
Is) = specific impulse (thrust /propellant flow rate) 
m = massratio = M)/M 
M = mass 
y = distance from center of planet 
t = time 
T = time parameter [Eq. (29)] 
V = velocity 
X = altitude parameter [Eq. (27)] 
Y = velocity parameter [Eq. (28)] 
a = gravity constant (distance*/time?) = 772g = roo*goo 
6 = flight path angle above local horizontal 
A = an increment or change in 
@ = range angle subtended at center of planet by vehicle’s 
motion 
— = approaches 
“A = equals by definition 
Subscripts 
c = circular value at initia 7 
0 = initial value 
00 = surface of earth 
M or min = minimum velocity 
E = parabolic or minimum escape velocity 


Received May 15, 1958. Revised and received September 4, 
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effective end of circular flight approximation 
due to thrust 


ll 


INTRODUCTION 


6 a SPECIFIC IMPULSE (ratio of thrust to propellant 
flow rate) of chemical rockets is decidely limited. 
Therefore, space flight systems employing chemical 
rockets will require large expenditures of propellant 
and, consequently, much bulky heavy hardware. On 
the other hand, propulsive systems developing ex- 
tremely high specific impulses such as ion rockets, solar 
powered vehicles,! and plasma jets will at least in their 
early development stages require large ratios of engine 
size or weight to thrust. This will limit the more effi- 
cient power plants to thrust accelerations much less 
than the value required to launch from the surface of 
the earth. Hence, their general use will be limited to 
phases of space flight commencing or terminating in 
satellite orbits. High-thrust chemical propulsive sys- 
tems will be used between the planet's surface and these 
satellite orbits. 

The purpose of this paper is to analyze and, insofar 
as possible, to generalize the low-thrust space trajec- 
tories in order to eliminate or at least reduce the amount 
of digital machine trajectory computation required for 
performance studies. 


(1) DescRIPTION OF TRAJECTORY 


When using a low tangentially directed thrust, the 
craft starting from a circular satellite orbit will follow 
an unwinding spiral trajectory, as illustrated in Fig. 1. 
Initially in the satellite orbit the gravity and opposing 
centrifugal acceleration are in balance. The application 
of a small amount of thrust along the flight path in- 
creases its centrifugal acceleration over that of the 
gravity acceleration. This causes the vehicle to move 
outward away from the planet, and the resulting con- 


Fic. 1. Symbols. 
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Fic. 2. Examples of oscillations in circular flight. 


version of kinetic energy to potential energy will re- 
duce the vehicle’s speed and consequent centrifugal 
acceleration below the opposing gravity acceleration. 
Therefore, a low-amplitude oscillation in velocity and 
altitude about the mean trajectory will exist. This 
oscillation damps out as the trajectory progresses. 

By definition the rate of curvature of the flight path, 
6, is the normal acceleration across the path divided by 
the vehicle’s velocity. Since the thrust is directed 
along the path, the only component of acceleration 
across the path is a gravity component. The effect of 
gravity diminishes as the vehicle moves away from its 
source, so the flight path ultimately approaches a 
straight line asymptotically. If this asymptote fails, 
as it does, to pass through the center of the planet, the 
offset is negligible compared to the distance of the 
vehicle at infinity, so the trajectory approaches vertical 
flight asymptotically. 

For these high specific impulse, low-thrust propulsion 


systems, the overall expenditure of propellant may be - 


of far less importance than the expenditure of food and 
supplies for the crew. For logistic reasons and to 
minimize the fatigue of the crew it may be desirable 
to reduce the overall flight time by maintaining thrust 
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far beyond the parabolic escape velocity. In some 
cases it may be advantageous to maintain power 
throughout flight, reversing it from an accelerating to 
a decelerating force in mid-course. 

Tangentially directed thrust was used throughout 
this study, because it results in the greatest instantane- 
ous rate of combined potential and kinetic energy in- 
crease. The analytical approaches presented here are 
equally applicable to the case of horizontally directed 
thrust and may even form the basis for subsequent 
studics of programed thrust direction. 


(2) CrrcuLAR FLIGHT OSCILLATIONS 
The general equation for radial acceleration is 


where the first two terms on the right-hand side are 
the centrifugal acceleration, the third term is the in- 
verse square gravity acceleration, and the last term 
is the vertical component of thrust acceleration. In 
the case of a low-thrust vehicle starting from a circular 
or near circular orbit, the vertical velocity, *, is com- 
paratively low, and the equation may be reduced to 
the following: 


= (V*/r) — (2) 


Where the overall changes are small, Eq. (2) may 
be put into the following incremental form by omitting 
higher order terms: 


At = (Vo?/ro) [1 + (2AV/Vo) — (Ar/r)] — 
(u/ro2) [1 — (2Ar/ro)] (3) 
where AV/Vp 1.0, Ar/ro < 1.0 (4) 


To eliminate the velocity variable from Eq. (3), the 
tangential acceleration, which is the sum of its thrust 
and gravity components, is used. 


V =f — (u/r’) (5) 
In incremental form for constant / this is 
AV = fAt — (u/ro?Vo)Ar (6) 


The equations are simplified somewhat by changing 
the independent variable from time, ¢, to the central 
angle, @. Within the validity of Eq. (4) we have 


(Vo/ro) At (7) 


Substitution of Eqs. (6) and (7) into Eq. (3) results in a 
linearized expression for radial acceleration. 


A? 4 d?Ar/dAg? = A + BAd — Ar (8) 
— (u/Vo), B= 2f(r0/Vo)? 
The solution of Eq. (8) is 
Ar = A(1 — cos Ad) + B(Ad — sin Ad) (9) 


The corresponding expression for velocity change 
is found by substituting Eqs. (9) and (7) into Eq. (6). 
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AV = (u ro? Vo) |} -A(1 — cos Ad) 
[1 — Vo2/2u)|BAd + Bsin Ad} (10) 


The constant, A, is seen from Eq. (8) to be the initial 
radial acceleration. Eqs. (9) and (10) show that this 
initial acceleration causes oscillations in altitude and 
velocity which go to zero at the end of each revolution 
about the planet. 

In the case of initial circular flight, the centrifugal 
acceleration balances the gravity acceleration. 


Au, ry” (11) 


The constant, A, then becomes zero, and further simpli- 
fication of Eqs. (9) and (10) leads to the following: 


Ar = 2f(ro/ Vo)? (Ag — sin Ag) (12) 
AV = —f(1/Vo) (Ad — 2 sin Ag) (13) 


Both of these equations show a sine wave oscillation 
about a mean path of constant slope. A graphical 
comparison of Eqs. (12) and (13) with a trajectory ob- 
tained by the stepwise integration of the classical equa- 
tions of motion on a high-speed digital computer is pre- 
sented in Fig. 2. It may be noted that, according to 
this simplified theory, the altitude levels off to hori- 
zontal flight at the end of each revolution and has its 
maximum rate of climb equal to twice that of the mean 
path at half a revolution about the planet. The ve- 
locity attains a maximum at one sixth of a revolution 
and a minimum at five sixths. Its steepest slope, oc- 
curring at one-half revolution, is three times that of the 
mean path while its slope at the start of each revolution 
is equal and opposite to the mean path. 

If the complete equation for radial acceleration, 
Eq. (1), instead of the reduced form, Eq. (2), had been 
used to derive Eqs. (9) and (10), and (12) and (13), the 
results would be much more complex equations of the 
same general form, but involving an exponential 
damping coefficient. The magnitude of these oscilla- 
tions actually damps out as the trajectory progresses. 
In the case of velocity, the peaks and dips become 
more equally spaced as their amplitude diminishes. 
It might appear from Fig. 2 that the amplitude of 
the altitude oscillation is increasing. This illusion is 
caused by the fact that the simplified theory gives a 
straight line for the mean path while in reality the 
mean altitude increases at an increasing rate. The 
time period of each oscillation of course increases be- 
cause the radius distance is increasing while the mean 
velocity is decreasing. 


(3) MEAN PaTH CIRCULAR FLIGHT 


The incremental equations for the mean path are 
obtained by omitting the trigonometric functions from 
Eqs. (12) and (13). 


Ar = 2f(r0/ Vo)? Ad (14) 
AV = —f(to/Vo)Ad (15) 


In accordance with Eq. (4), the smaller these incre- 


ments become the more valid the equations. We now 
allow these increments to approach zero and become 
derivatives. The constants denoted by subscripts 
zero then become the local value of the variables. 


dr = 2f(r/V)? do (16) 
dV = —f(r/V)do (17) 
Dividing Eq. (16) by Eq. (17) yields the differential 
equation of the mean path, 
dr/dV = —2(r/V) (18) 
which is readily integrated to the following: 
(r/o) (V/Vo)? = 1 (19) 
Since the vehicle started in circular flight, we have 
from Eq. (11), 
roVo? = (20) 
Therefore Eq. (19) shows that the approximate mean 
path remains in the circular velocity relationship. 
Eventually, of course, gravity will diminish to the point 
that the thrust acceleration will destroy this equilibrium 
and the spacecraft will eventually accelerate without 
limit. 
Substitution of the differential form of Eq. (7) into 
(17) and integration yields 


V, Vo = 1 = (Vr; Vo) (21) 
where the ideal or thrust velocity is 
t 
Ve sat (22) 
0 


Eq. (21) indicates that the decrease in vehicle velocity 
is numerically equal to the velocity contributed by 
thrust. This means that the velocity lost to gravity is 
exactly twice that gained due to thrust. 

The equation for distance from the center of the 
earth is found by substitution of Eq. (21) into (19). 


= {1 — (Vr/ Vo) (23) 


Eqs. (21) and (23) express velocity and radius dis- 
tance in terms of the thrust velocity, Eq. (22). In the 
special case of constant specific impulse with either a 
constant or variable thrust acceleration /, the thrust 
velocity may be given in terms of the mass ratio. 


Ver = In m (24) 


In that case, the following simple relation holds be- 
tween time and mass ratio: 


t = (Lspgo0/fo) [1 — (1/m)] (25) 
In any case, where f is constant, 
Vr = ft (26) 
(4) GENERALIZED TRAJECTORIES FOR CONSTANT 
THRUST ACCELERATION 


There are at least two cases of practical importance 
where constant thrust acceleration may be assumed for 
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Fic. 3. Parametric velocity vs. radius distance. 


trajectory computation. One example is a low-thrust 
spacecraft with a specific impulse so high that the over- 
all expenditure of propellant is small compared to the 
vehicle’s weight as it leaves the satellite orbit. An- 
other example is a higher-thrust vehicle operating con- 
tinuously at a maximum allowable thrust acceleration. 
Such a vehicle would require throttling of the thrust as 
the overall mass diminished due to propellant expendi- 
ture. The first step in this analysis is the generalization 
of the equations of motion. 


(4.1) Generalized Equations of Motion for Constant 

Thrust Acceleration 

Radial and tangential acceleration are given in Eqs. 
(1) and (5), respectively, as functions of V’, 7, ¢, u, and 
f. These basic differential equations may be expressed 
in a form independent of both the gravity constant and 
the thrust acceleration by the introduction of three 
dimensionless parameters. The derivation of these 
parameters is not presented here, though their validity 
is self-evident. The “altitude,” ‘‘velocity,”’ and ‘‘time’’ 
parameters, X, Y, and 7, respectively, are defined as 
follows: 


X = a*(r/ro) = (f/m)? 1 (27) 

¥ (V/V) = (uf) (28) 

T = (Ve/na t= (29) 

where a = f/(u/r?) (30) 


The dimensionless thrust acceleration, a, is in units of 
gravity (g) at the radius, ro, of the initial satellite orbit. 
In the special case of trajectories starting from circular 
orbit, V, = Vo (see Symbols). 

Substitution of Eqs. (27), (28), and (29) into Eqs. (1) 
and (5) results in the generalized differential equation 
for radial and tangential acceleration in terms of X, Y 
and their derivatives with respect to the dimensionless 
time parameter, 7. 


X = (¥?/X) — (X?/X) — (1/X2) + (X¥/Y) (31) 
(32) 


These equations may easily be combined to com- 
pletely eliminate the time parameter; resulting in an 
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equation of X and Y and the first and second derivative 
of Y with respect to X. 

It must be emphasized that these are differential 
equations and the final curves resulting from their in- 
tegration will in general depend upon the initial values 
of X and ¥ and the first derivative of XY with respect to 
Tor Y. Inany instance where the curve is independent 
of the initial conditions, one curve will suffice for all. 
Such an instance occurs when the circular flight de- 
scribed in Section (3) and Eq. (19) exists. In that 
event, the vehicle will have the same mean velocity, 
altitude, and flight path angle as it passes a given 
point, regardless of where or when the initial point of 
the trajectory may have occurred. In this case, one 
curve of X versus Y or X or Y versus 7 may be ex- 
pected to apply for all trajectories on out to infinity. 


(4.2) Asymptotes 


The approximation for the mean circular path, Eq. 
(19), may be put in terms of Y and_X as follows. 
y bey 


= 1 (33) 


In terms of logarithms of Y versus X this is a straight 
line with a slope of —1/2 passing through the origin. 
The derivation of this equation stemmed from the 
assumption that the vertical velocity, *, was negligible 
and was hence set equal to zero in Eq. (1) to obtain 
Eq. (2). Actually the mean path must have some ver- 


tical velocity or the mean path altitude would not 


change. The sine of the flight path angle is the radial 
velocity over the tangential velocity and may be found 
by combining the differential of Eq. (33) with Eq. (32). 


sind 4 ¥/Y = 2X? = 2/Y! (34) 


The only place this mean angle would equal zero pre- 
cisely is at X = OorlogX = —o. Therefore Eq. (33) 
represents an asymptote. Since Eq. (34) was obtained 
from Eq. (33), it too represents a straight line asymp- 
tote for sin 6 when plotted on double log scale. 

The slope of the circular asymptote for a plot of Y vs. 
T is obtained by differentiating Eq. (21) with respect 
to time, /, and putting the result in dimensionless form. 


dY/dT = —1 (35) 


During the circular phase of flight, the velocity 
diminishes approximately linearly with time. Even- 
tually gravity diminishes so much that the mean path 
departs appreciably from its circular asymptotes and 
a minimum velocity will be reached beyond which the 
vehicle continues to gain in speed without limit. If 
thrust is maintained long enough, the gravity will be- 
come negligible and, as explained in Section (1), the 
trajectory will approach the vertical asymptotically. 
Under these conditions we might expect the velocity, 
altitude, and time relationships to become asymptotic 
also. The equations for the flight path at infinity will 
also be the equations of these asymptotes. 

Because, as the vehicle approaches infinity, gravity 
approaches zero, the flight path angle approaches the 
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vertical, and time approaches infinity, the following 
relations must be true: 


(36) 
vos (37) 
r—> (1/2)/? (38) 
V— ft (39) 


From Eqs. (27) and (28), 


din Y r 


(40) 
din X Vi 


Substitution of Eqs. (36) through (39) into Eq. (40) 
yields the logarithmic slope of the X — Y infinity 
asymptote. 


din Y 1 


= (41) 
din X 2 


Since the slope of the corresponding asymptote for 
circular flight was minus 1,2, this might suggest that 
the mean path of Y vs. X is a hyperbola in logarithms. 
Such is not the case, though it may be possible that the 
mean path is the sum of such a hyperbola and a long 
period oscillation. This has not as yet been verified, 
nor has a simple algebraic expression for the entire 
mean path been obtained. 

The final equation for the infinity asymptote is 
obtained by substituting Eqs. (36) and (37) into Eq. 
(40), putting the results in terms of Y and Y and 
equating them to Eq. (41). 


=2 (42) 


Unlike the circular asymptote, Eq. (33), validity 
of this infinity asymptote, Eq. (42), is not limited to 
vehicles using low thrust and starting in circular flight. 
All trajectories involving constant tangentially directed 
thrust acceleration would tend to eventually approach 
Eq. (42). 

The slope of the velocity-time infinity asymptote is 
found by putting Eq. (37) in dimensionless form. 


dV/dT = 1.0 (43) 


(4.3) Generalized Trajectory Plots 


Fig. 3 is a plot of several typical trajectories in para- 
metric form of Y vs. X. Straight dashed lines indicate 
the circular and infinity asymptotes and the locus of 
minimum escape velocity points. The initial points 
of the trajectories are circled and are easily located as 
occurring at X = a‘ because r/7p is initially unity. 

These trajectories were obtained in the traditional 
manner from a high-speed digital computing machine 
and then plotted in the dimensionless parametric form 

It should be noted that all trajectories starting from 
circular orbits (along circular asymptote), with an a 
< 10~°, fall on the same curve. This single curve has 
a minimum velocity point and a parabolic velocity point 
at definite values of X and Y. It crosses the infinity 
asymptote at an XY of about 1.5 and then proceeds to 
approach it from below. As shown, higher acceler- 
ation trajectories such as 10~' and 10° g starting in 
circular flight overshoot this mean curve and then 
eventually fair into it. Still higher accelerations such 
as an a of 10 or 100 initial g simply fair into the mean 
curve. 
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Mean path parametric velocity vs. radius distance. 
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Fic. 5. Velocity parameter vs. time parameter from 
escape point. 


Trajectories of low acceleration (a < 10~*) starting 
with velocities a few hundred feet per second above 
circular velocity will oscillate about the mean path and 
rapidly merge with it (examples not shown). 

One example of a trajectory with an acceleration, 
a, of 10~* and an initial velocity of twice circular (well 
above escape) is shown in Fig. 3 to approach the in- 
finity asymptote. 

Fig. 3 shows that for the higher accelerations, say 
10*+!, the curve is almost vertical between circular 
and parabolic velocity. This means that for purposes 
of estimating time to reach minimum escape, the 
change in altitude may be ignored and the time is then 
simply the difference between escape and circular 
velocity at the initial altitude divided by the thrust 
acceleration, f. In terms of our dimensionless param- 
eters, this is 


Ty = a4 1) (44) 


A more precise approximation for Eq. (44) is pre- 
sented in reference 2, though the difference is less than 
one part in two hundred for values of a greater than 
one. 

For computational purposes the curve for low acceler- 
ation initial circular flight trajectories from Fig. 3 is 
reproduced to a larger scale in Fig. 4. A corresponding 
single-valued curve should exist for the case of hori- 
zontally directed low thrust, starting from circular 
orbits. It may be presented as part of future analysis. 
In reference 3, Tsien dealt with the two cases of hori- 
zontally and radially directed thrust, but not tangen- 
tial thrust. His work was confined to determining 
conditions at only the parabolic escape point. His re- 
sults at this point for the case of horizontal thrust are 
not greatly different from those presented here for tan- 
gential thrust. 

Fig. 5 is a parametric plot of mean velocity vs. time 
for the low constant tangential thrust acceleration 
trajectories of Fig. 4. Time is presented as the differ- 
ence between flight time and time required to reach 
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“TABLE 1 
Vi Vin Ve 
¥ 1.85 1.494 1.509 
6 9.4° 39.2° 
T a~/4.1 85 960 1-1/4. 809 


escape velocity. If, instead of escape velocity, the 
minimum velocity or the time to achieve any arbi- 
trarily fixed value of X had been selected as the sub- 
tractive quantity, the same curve would exist with 
only the abscissa scale being shifted. The asymptotic 
slopes of this plot are #1 in accordance with Eggs. (35) 
and (43). Initial points on this curve are of course 
located at Y = a~'/*. Curves of Y vs. T — Tex for 
accelerations greater than 10~° g for trajectories start- 
ing from circular orbits are presented in Fig. 5 as 
dashed lines. 

It may be shown that the number of revolutions 
minus the number to any fixed point such as the para- 
bolic escape point is a single-valued function of any of 
the coordinates of Figs. 4, 5, 6, or 7. 

Parametric velocity and radius distance as functions 
of the sine of the mean flight path angle for these very 
low acceleration trajectories are presented in Fig. 6. 
The circular flight asymptotes of Eq. (84) and the in- 
finity asymptotes are also shown. The corresponding 
plots of X and Y vs. 6 in degrees are presented in Fig. 
7. It is interesting that minimum velocity and para- 
bolic velocity occur at fixed flight path angles regard- 
less of the values of /, uw and 7%. 


(4.4) Specific Solutions 


Specific solutions for the velocity, altitude, flight 
path angle, and flight time to any point on the mean 
path of any trajectory using constant low (a < 10~*) 
tangentially directed thrust and starting from a circular 
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orbit may be obtained from Figs. 4, 5, and 7. Specific 
solutions for the effective end of the circular flight 
approximation of Section (3) and for the minimum 
velocity and minimum escape velocity are presented 
in Table 1. 

Using the following values of earth radius and 
gravity constant, and expressing f in units of ft./sec.’, 
yin feet, V in ft./sec., and ¢t in seconds, some of these 
equations for the planet earth are presented below. 


= 2:08555 10° ft. 


At minimum velocity: 


ru = 8.78 X 107 f-!/? 

Var = 1.626 X 10* f*/* 

64 = 33.3° 

ty = (Vo = 1.045 104f1/4) 
At parabolic velocity: 

re = 1.042 X 

Ve 1.643 104/'/4 

39.2° 

te = f-'(Vo — 8.81 X 10%f!/4) 


A purely analytical approximation of time to reach 
parabolic velocity with very low thrust is presented in 
reference 2. As compared with the accurate digital 
machine results contained here the method of reference 
2 yields a time that is 7 per cent low for an acceleration, 
a, of 10~4 and is 28 per cent low for an a of 10~*. This 
is very good for a purely analytical approach, but, 
unfortunately, the method used cannot yield any 
information about the trajectory other than the time 
to escape. 


(5) CONCLUSIONS AND RESULTS 


The trajectory for extremely low-thrust vehicles 
starting in circular orbit consists, for some time, of a 
damped oscillation about a mean path. The circular 
velocity-altitude relationship holds for this mean 
path. The mean velocity is equal to the initial velocity 
minus the ideal or thrust velocity. Since the thrust 
velocity is the integral of thrust acceleration with 
respect to time, the equations for velocity and distance 
from the center of the earth may be put in terms of time 
or mass ratio. These results are summarized in Eqs. 
(19) through (26). 

In the case of constant thrust acceleration circular 
flight, the approximate velocity and altitude oscillations 
are presented in Eqs. (12), (13), and (7). 

By use of the parameters of Eqs. (27), (28), and (29) 
the differential equations of motion may be generalized 
for all cases of constant thrust acceleration regardless 
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of thrust direction. This generalization makes the 
differential equations independent of the thrust acceler- 
ation and the gravity constant (choice of planet). It 
has applications beyond the scope of this paper. 

Applying these generalized parameters to the special 
case of tangential constant thrust acceleration from 
circular orbits results, for the case of thrust acceleration 
of 10~' or greater initial g, in a separate curve relating 
the velocity parameter to the altitude parameter or 
time parameter for each acceleration. These curves, 
Figs. 3 and 5, apply for any initial altitude and the 
gravity constant. For those trajectories of this group 
using thrust accelerations of 10~* or less initial g, an 
even greater generalization is possible and one curve 
results for all thrust accelerations (Figs. 3, 4, 5, and 7). 
As shown in the Figures and Eqs. (33), (35), (42), and 
(43), these curves approach very simple known asymp- 
totes in their early and late stages. 

Unique solutions such as those for the parabolic 
velocity point are apparent (see Section (4.4) and 
Fig. 4). 
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The Problem of Obtaining High Lift-Drag 


Ratios at Supersonic Speeds 
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INTRODUCTION 


HE IMPORTANCE of the lift to drag ratio is well 

known to all aircraft designers since it gives, to a 
great extent, the aerodynamic efficiency of the air- 
plane. Aerodynamic efficiency, however, is only one 
component of the grand compromise that a completed 
airplane represents. At subsonic speeds, lift-drag 
ratios of well over 200 have been measured in wind 
tunnels on airfoil sections; but few powered aircraft 
have attained (L/D) value of 20. It is invariably true 
that the requirements of stability and control, struc- 
ture, and flight operation all contribute to reducing the 
design (L/D)ma, considerably below those exotic values 
which can be predicted from unrestricted aerodynamic 
theory. If, however, a certain range or operating 
efficiency is required, there is most certainly a minimum 
(L/D)maz Value for which the goals are just attainable. 
If we examine the range equation we see that range is 
proportional to the lift-drag ratio, the thermopropulsive 
efficiency, and the logarithm of the initial to final 
weight ratio. The appearance of the lift-drag ratio as 
a linear factor in the range equation indicates that every 
attempt should be made to increase (L/D) maz; how- 
ever, the search for higher (L/D),.,; may lead to strange 
and unorthodox configurations. Most frequently, such 
configurations are ruled out by the adverse effects of 
their geometry on the weight ratios. In the present 
paper, we will deal with the maximum lift-drag ratio 
problem for conventional configurations having a wing 
and a body in close proximity to each other. No at- 
tempt will be made to select a particular configuration 
as being the best. However, the promising direction 
to go from the aerodynamic view will be stressed with 
the understanding that the other factors may outweight 
the aerodynamics. 


ANALYSIS 


General Equations 


The assumptions of the analysis are as follows: (1) 
that the skin-friction drag as well as the wave drag due 
to thickness is independent of angle of attack, (2) that 
linearized theory is adequate for performance estimates, 
and (3) that only a small error arises from the simple 
addition of component wave drag due to thickness. 
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The last assumption is probably valid at Mach Num- 
bers not too close to unity. The theoretical drag polar 
for any symmetrical configuration may then be ex- 
pressed as follows: 


Cp = Cp, + C1? (1) 


where Cp, is the drag coefficient at zero lift and is 
composed of both wave and friction drag, (OCp/OC;”) 
is the drag rise factor which is theoretically a constant, 
and C, is the lift coefficient. For nonsymmetrical (with 
respect to the lift direction) configurations, the polar is 
still a parabola but is, on the basis of theory, displaced 
with respect to the origin. When a group of configur- 
ations having a common form but varying degrees of 
asymmetry are considered, a so-called ‘“‘rubber airplane” 
drag polar may be obtained in which lift coefficient rep- 
resents the design state and the polar is symmetric 
about the vertical axis. 

The solution of Eq. (1) for maximum (L/D) gives 
the following simple relation: 


(L/D) maz = (1/2) 


which is valid for either the symmetric or the “rubber” 
asymmetric configuration. Attainment of high 
(L/D) max Values is thus seen to depend on the two fac- 
tors Cp, and (OCp/OC,”). The drag rise factor (OCp + 
OC,”) is primarily affected by the wing selection, 
whereas the term Cp, is influenced by several factors— 
skin-friction coefficient, wave drag due to volume, and 
ratio of total wetted area to wing area. The average 
skin-friction coefficient may be greatly reduced by ob- 
taining a substantial run of laminar flow. However, 
experience has proved that laminar flow is extremely 
difficult to retain at the high flight Reynolds Numbers 
without extremely smooth surface construction and 
continued meticulous maintenance as well as boundary- 
layer suction. The attainment of laminar flow at high 
Reynolds Number on operational aircraft is certainly 
to be desired, but has so far been an unfulfilled obsession 
of the aerodynamist. We have therefore resigned our- 
selves in the present analysis to completely turbulent 
flow for which the friction coefficient is determined by 
the Reynolds and Mach Numbers. 

With fixed friction coefficient, the fuselage or fuselages 
should be as close to spherical as the wave drag will 
allow. Of course, the wave drag at supersonic speeds is 
intense, and, hence, our tendency toward spherical 
shape is persuaded to settle for a fineness ratio of order 
1/15. It is important to note that two optimum 
bodies of revolution having no mutual interference 
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effects produce about 25 per cent more drag than a 


single optimum body of the same total volume. The 
point we are leading up to here is that, for given volume 
requirements, the drag of the fuselage is pretty well 
fixed, so that the choice of a wing becomes the main 
problem. The wing not only controls the drag rise 
factor (OCp/OC,”) but has inputs in Cp, by virtue of 
its own wave drag and its influence on horizontal and 
vertical control surface size. It is true though that 
the effects on Cp, are secondary with respect to the 
drag rise factor since the variations in, say wing drag, 
are only a friction of the total which then appears 
under a square root sign in the (L/D), relation. 

Let us turn then to the question of how the drag due 
to lift can be reduced. It is true that, for uncambered 
wings having sweepback such that the flow components 
normal to the leading edge are greater than sonic 
(supersonic edges), the drag rise factor is inversely pro- 
portional to the lift curve slope. Now the lift curve 
slope of a two-dimensional (infinite aspect ratio) wing 
is 4/V M2? — 1 and so is it also for all supersonic-edged 
delta wings. The only way the lift curve slope can be 
raised is by carving out the rear center of a delta 
wing, creating an arrow wing and thus increasing the 
aspect ratio. The greatest gains occur when the lead- 
ing-edge sweep necr3 the sonic edge condition. Such 
high aspect ratio wings are not too attractive, however, 
because the wave drag tends to become excessive for 
realistic structural weights. For the supersonic lead- 
ing-edge condition, there are still some other ways for 
reducing the drag due to lift. 

The first method is to employ a camber (wing warp) 
which provides a favorable load distribution. This 
technique has been developed by the methods of linear- 
ized theory by Jones, Graham, Ward, and others! 
and has led to the idea of a ‘‘lower bound”’ in drag due 
to lift for a wing of given plan form. The work of 
Zhilin® and Germain’ has shown that the gains to be 
obtained over the drag of a flat wing are maximum at 
the sonic edge condition and further that the optimum 
sonic-edged delta wing theoretically produces a drag 
due to lift only 11 per cent lower than the flat wing 
of similar plan form. In terms of (L/D), then a 
6 per cent increase in (L/D); may occur from the 
camber. 

Additional gains may be made by carrying lift on 
the fuselage (described in reference 8). In this manner 
the loading is distributed fore and aft such that the 
wave drag component of the drag due to lift is reduced. 
The amount of lift to be carried on the fuselage is lim- 
ited by the effects of separation and shock formation, 
so that it is difficult without a careful experimental 
study to determine the gains to be realized. It should 
at very least be possible to offset the trim drag by carry- 
ing lift on the fuselage. 

A second method for reducing drag to lift of super- 
sonic edged wings which has received much attention 
is the so-called interference method, whereby the wing 
is situated in the pressure and upwash field of a body.*~"! 
In much of the published work, it is hard to assess the 
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value of the interference effect because of the difficulty 
in finding a proper configuration for comparison. We 
have therefore studied a simplified wing-body inter- 
action configuration from which a comparative study 
is easily made and which included the effects of base 
pressure and skin friction. The model studied is 
shown in Fig. 1. It consists of a sonic-edged delta wing 
of zero thickness and a parabolic-are body of revolution, 
or half-body, as the case may be. For the interference 
study the half-body is placed under the wing as shown 
and two parameters are varied in an attempt to obtain 
a fairly wide variation in the induced pressure distribu- 
tion under the wing. The fineness ratio of the basic 
shape was varied as well as the base cutoff point, ex- 
pressed as a fraction of the basic shape length, /. 
Thus, the body shape varied from a full parabola to a 
cone. Curves are shown for the cases in which the 
volume was held constant for a given wing area as 
shown by the curves of constant V*°/S. A definite 
maximum in (L/D) n,; is found for each volume ratio. 
Note that each point represents a different body shape. 
There is a marked effect of the base drag in these results 
and it may be argued that some of the base area can be 
filled by the expanded jet exhaust. This is true for 
nestled engine packs; however, the duct also removes 
usable volume, and the computation becomes rather 
involved. We have, however, computed the effects 
for no base drag at all at JJ = 3, and a summary of the 
results is shown in Fig. 2. A comparison is made 
here for bodies of similar volume and base area; the 
Figure is nearly self-explanatory except that the curve 
labeled symmetrical body-cambered wing represents the 
case of a symmetrically disposed body on an optimum 
cambered wing, thus claiming the 6 per cent (L/D) maz 
increase mentioned previously. Note that the (L/D) maz 
values are unrealistic because no drag due to wing 
volume, jet engine installation, control surfaces, and 
trim have been included. The lesson to be learned 
here is that the magnitude of the favorable interference 
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Fic. 6. Oil film flow picture of arrow wing; fixed transition, C, = 0.1, W = 2.87. 
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effect is of the same order as that to be obtained by 
wing warp, regardless of the base drag. We see then 
that 6 per cent gain in (L/D), is of the order of the 
increase to be anticipated from sonic-edged wings by 
the use of camber or interference effects. The authors 
feel that the gains theoretically obtainable by the use 
of interference or camber are not additive, since in each 
case the gain would appear to arise from the attain- 
ment of a more nearly ideal load distribution. 

The linear theory indicates that rather large reduc- 
tions in drag due to lift as well as wave drag due to 
thickness should be obtained by use of sweptback wings 
with subsonic leading edges. This result has been 
known for some time due to the research of R. T. 
Jones'* and others.'*:'* For flat wings the theory 
predicts a favorable leading-edge thrust which unfor- 
tunately is rarely found to any appreciable extent in 
the experimental results. For this reason the flat 
swept wings lose much of the predicted drag-reduction. 
In cases where the trouble is associated with leading- 
edge separation, tests have shown that conical camber 
or leading-edge camber has allowed some of the theo- 
retical performance to be recaptured.'*: Recently, 
however, work by various investigators previously 
mentioned has shown that, by optimizing the wing load 
distribution, a camber surface is obtained which requires 
no leading-edge suction at all and yields theoretical 
drag due to lift values slightly below those of the flat 
wing with full leading-edge suction. The results of 
such theorizing are shown in Fig. 3. Here the cam- 
bered wing values represent results obtained by super- 
position and optimization of 5 lift loadings. The 
improvement over the flat wings is small, but the im- 
portant thing to remember is that no leading-edge 
suction is demanded. 

Various techniques for calculating the optimum have 
recently appeared.®: *:'” It will be noted that these 
methods have in common the assumptions of linear 
theory; they differ considerably in method but should 
obtain the same drag values and wing camber surfaces. 
The Figure clearly indicates the large gains theoretically 
available by the use of high aspect ratio wings with 
subsonic leading edges. For example, the arrow wing 
with a 40 per cent cutout at a value of 8 cot A = 0.5 
theoretically provides a 39 per cent reduction in drag 
due to lift over a flat sonic-edged delta wing. 

We come now to the question of the theory’s validity. 
The linear theory applied to swept wings has in it the 
assumption of constant Mach angle. That is, all 
disturbances can propagate as far forward as the fixed 
Mach cone originating at the disturbance. In reality 
the local Mach lines may vary widely in their inclina- 
tion to the free-stream direction. The terms subsonic 
and supersonic edges derive from the idea of fixed Mach 
angles, and the nature of flow development around a 
leading edge depends uniquely on this assumption. In 
the real flow the character of the local Mach lines is 
determined by the pressure and flow direction such 
that on a so-called subsonic edge the local Mach lines 
may lie behind the edge. In this case, the local flow 


component normal to the leading-edge line is larger 
than the local sonic speed and the wing is effectively in 
a transonic flow field even though the free-stream 
Mach Number is far above one. Although the non- 
linear equations of motion cannot yet be handled, it 
is clear that when the local flow becomes ‘‘transonic’’ 
the pressure development on a wing will deviate mark- 
edly from that predicted from linear theory. Since 
an ‘“optimum’’ wing camber surface requires a careful 
balancing of wing slope and pressure, the deviation due 
to transonic effects will certainly cause a rapid drag 
rise similar to that found experimentally on two- 
dimensional cambered airfoil sections as the critical 
speed was exceeded. It is possible to make a crude 
calculation of the onset of ‘‘transonic’’ flow over swept- 
wing leading edges by using the local pressure coeffi- 
cient on the wing upper surface and the local flow in- 
clination computed from simple sweepback theory. 
We then find a critical Mach Number, which is a func- 
tion of the sweep angle, the pressure loading at the 
leading edge, and the airfoil section. A plot of this 
relationship is shown in Fig. 4 for the case of uniform 
load. Optimum loadings or flat wing loadings increase 
the negative pressure on the leading edge and hence 
tend to reduce the critical Mach Numbers from those 
shown. 

The Figure shows that rather large sweep angles are 
required at Mach 2 or above if a subsonic type of flow 
is to be maintained over the wings. We note that 
at M = 3 approximately 80° of sweep is required just 
for the case of uniform load at design C, of 0.1. There 
is an important effect of leading-edge profile, and it 
would appear that gains may be expected from the use 
of sharp-edged airfoils. A recent test has shown that 
a severe penalty may be incurred for exceeding the 
critical Mach Number. Fig. 5 shows a cambered 
arrow-type wing having 75° of sweep and designed for 
nearly optimum drag due to lift at 1/J = 3. The theo- 
retical drag rise factor is 0.44, but in test it produced 
over 0.7. The high drag is most probably caused by 
the large deviation of the actual upper surface pressure 
distribution from that of the predicted one which was 
a linear variation. The small dashed lines on the pres- 
sure plots show the pressure coefficient at which sonic 
normal velocity component is reached. Clearly the 
flow is well into the transonic range. The overexpan- 
sion produces shock waves on the upper surface which 
are inclined to the free stream at angles less than the 
leading-edge angle. The shock system produces a 
boundary-layer separation in the cross flow which can 
be seen in Fig. 6. This photograph, taken in the Lang- 
ley Unitary Plan Wind Tunnel, shows the half wing 
in the tunnel with flow from left to right. The flow 
is attached to the upper surface up to the first white 
line; thereafter the flow separates and rolls up into a 
flat vortex. 

Unfortunately, we cannot show any conclusive 
evidence that cambered wings which are kept below 
their critical speed produce the optimum drag rise 
factors predicted by theory. However, tests of a 68° 
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sweptback delta wing cambered for uniform load at 
M = 1.62 (see reference 13), have indicated a laminar 
flow on the upper surface free from separation from lead- 
ing to trailing edge. The flow over the wing was judged 
to be well below critical conditions. 

The conclusion to be drawn from this discussion is 
that sweptback wings designed according to linearized 
theory for optimum drag due to lift must be carefully 
checked to be sure the design is not too close to critical 
conditions of flow. It is also necessary to check the 
design pressure gradients in the cross-flow direction 
and the magnitude of the trailing-edge pressure dis- 
continuity to be sure that the boundary layer is not 
too close to separation. These checks on the boundary- 
layer separation are difficult to tie down because there 
are little or no data available on turbulent flow pressure 
characteristics in the presence of a strong lateral flow. 
Nevertheless, simple sweep or cross-flow analysis 
should provide the proper order of magnitude. 


CONCLUDING REMARKS 


Returning to the main idea of obtaining high 
(L/D) mar, it may be instructive to put some simple 
estimates into the (L/D), equation. If we assume 
a M = 3 airplane having a wetted area to wing area 
ratio of 3, a turbulent friction coefficient of 0.0014 
corresponding to a high Reynolds Number, we obtain 
approximately 0.0040 for the friction drag coefficient. 
A wing and tail surface wave drag of about 0.0015 
would be realistic. A fuselage wave drag of 0.0010 
and air inlet and exit drag of 0.0010 brings us to the 
“ball park’’ figure of 0.0075 for Cp,. Use of a super- 
sonic-edged delta wing would give a drag rise factor of 
0.7. Thus an (L/D), value slightly under 7.0 is ob- 
tained at 1/ = 3 assuming no penalty for trim. Most 
designers will see at once the great difficulty in reducing 
the Cp, by a sensible amount since Cp, has a typical 
tendency for growth once the design has left the paper 
stage. For this reason then, the drag rise factor is 
emphasized as being the best hope for a substantial in- 
crease in (L/D),.,,; above the quoted values. Further 
manipulation of the numbers in the (Z/D),,,, equation 
is left to the reader. 

To summarize the foregoing discussion, we have 
attempted to show that high aspect ratio sweptback 
wings still offer the best hope for substantial reductions 
in drag due to lift. Sweptback wings which are cam- 
bered and twisted for optimum loading have been 
shown theoretically to provide, without leading-edge 
suction, a drag due to lift slightly lower than flat wings 
with full leading-edge suction. Finally some real 
flow limitations to the pressure loadings and sweepback 


angles have been pointed out and a new supersonic 
range critical speed has been stressed. 
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Landing-Gear Vibration Instability 


FREDERICK L. RYDER* ann MELVIN ZAID* 
Republic Aviation Corporation 


SUMMARY 


A theory of gear walking, or excessive vibration of the landing 
gear during the landing run, is given which attempts to isolate 
the essential physical parameters and to derive the differential 
equations which must be satisfied by the resulting idealized model 
of the aircraft. It is found that if the curve of wheel-braking 
torque vs. angular velocity has a negative slope under conditions 
of rapid variation of the angular velocity (as might be expected 
on the basis of conventional friction data), then the model will be 
unstable below a certain critical aircraft speed. While this speed 
is, in general, relatively low, in high-speed aircraft with relatively 
flexible landing gear and connections it may be high enough for 
the landing gear vibration, resulting from the instability, to reach 
objectionable levels with consequent gear walking. 

If, on the other hand, the slope of the curve is zero or positive, 
then the model has a definite margin of stability, so that gear 
walking can occur only if the braking torque repeatedly exceeds 
the skidding value. At low aircraft speed, where the margin of 
braking stability becomes small, overshoot of the tangential 
force between wheel and ground due to initial braking disturb- 
ance may aid in causing repeated skidding. 

A detailed study is made of the behavior under sinusoidal dis- 
turbances for the purpose of gaining a physical understanding 
of the phenomenon. Remedial measures are discussed briefly, 
including the application of damping devices and the manner 
in which the critical speed may be lowered by varying certain 
parameters. Experiments for obtaining data required in the 
practical application of the present theory are suggested. 


INTRODUCTION 


ct VIBRATION of the landing gear during the 
landing run, or gear walking, has been observed in 
aircraft where the gear and its structural connection 
with the fuselage are unusually flexible. The phenome- 
non, which is usually observed toward the end of the 
landing run when the aircraft speed is relatively small, 
can be sufficiently severe to set the entire fuselage into 
violent oscillation, and has been known to cause 
structural failure. 

Wignot and Hoblit! have analyzed the conditions for 
landing gear instability with the wheel assumed 
locked. Their results do not shed light directly on the 
usual phenomenon in which the braking torque is not 
large enough to cause wheel locking. 

Others have suggested that a periodic forcing func- 
tion results from slight irregularities of the braking 
surfaces, since at certain low speeds these irregularities 
can be shown to lead to a periodic disturbance of the 
same frequency as that observed in gear walking. 
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Similarly, certain landing field irregularities can be 
found which will introduce periodic disturbances at 
the gear walking frequency at the speeds of interest. 
However, it is considered here that such sources of dis- 
turbance are too uncertain to be of significance in con- 
nection with gear walking, and, instead, we shall seek 
an explanation in a more detailed study of the me- 
chanics of braking of an aircraft with flexible landing 
gear. For this purpose, we shall attempt to make those 
assumptions and idealizations which reduce the phe- 
nomenon to its essentials, so that, even though our re- 
sults may not be quantitatively accurate, they may at 
least give an adequate qualitative understanding of 
gear walking. 


SIMPLIFIED MODEL OF AIRCRAFT 


The aircraft is considered perfectly rigid except for 
the landing gear and its structural connection to the 
fuselage, which are assumed to be flexible only as re- 
gards fore and aft movements of the center of the 
landing wheel, this flexibility being indicated by the 
spring of stiffness k in Fig. 1. Let M indicate a fictitious 
mass which, when placed at the wheel center, would 
result (with actual masses neglected) in the same period 
of free vibration as that of the fundamental mode of 
fore and aft vibration of the wheel center with actual 
masses present. Thus, J/ includes the mass of the 
wheel and part of the mass of the landing gear and its 
connections. 

The fuselage is considered to have infinite mass and 
to be traveling at velocity |’; for simplicity we con- 
sider the fuselage fixed in space and consider the ground 
to be traveling at velocity V’, as indicated in the sketch. 
Fig. 1 also shows the brake torque 7, the ground force 
F, the wheel angular velocity w, and the translational 
displacement x of the wheel center, the latter measured 
from the position corresponding to zero force in the 
spring. The normal force between wheel and ground, 
which is assumed constant, is not shown. 

The effects which have been omitted from Fig. 1 
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Fic. 1. Simplified model of aircraft. 
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Fic. 2. Ground force F vs. slip S. 
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include such things as the vertical spring compliance 
of the tire and landing gear, tire hysteresis and cant 
angle, and the fact that the wheels on both sides of the 
aircraft do not act in perfect unison. While such 
effects may lead to important practical consequences, it 
is felt that their inclusion here would tend to mask the 
phenomenon of present basic interest. 

Force F is a function of w, V, the elastic and fric- 
tional characteristics of wheel and ground, and the 
normal force between wheel and ground. Experi- 
mental data given in Figs. 53 and 54 of Smiley and 
Horne’ show that, as w and V alone vary, F is deter- 
mined by the slip (for steady values of w and V’), de- 
fined thus in the absence of spring deflection: 


S = (V — Rw)/V (1) 


where R, is the effective wheel radius—that is, Rw = 
V—in free rolling (fF = 0). For a given wheel and 
landing field, and a prescribed normal load, the relation- 
ship between F and S has the general form shown in 
Fig. 2, with a peak at S,, and F,,. An intuitive appre- 
ciation of this curve can be gained by considering Eq. 
(1) under the following conditions: 

(1) If the wheels are free and F = 0, then S = 0 
because Rw = V. 

(2) If the wheels are locked and w = 0, then S = 1, 
and F is appreciably lower than its peak value. 

(3) As ground force F rises from zero, it acts in con- 
junction with the equal and opposite shaft force to 
stretch the tire so that a single revolution of the shaft 
results in a greater ground distance traveled; this cor- 
responds to a decrease in w, or a rise in S. Hence, S 
continues to rise with increasing F until the latter 
reaches its peak, which occurs when 5S is of the order 
of 0.2; as S increases, slippage commences between the 
wheel and ground, and becomes greater and greater, 
leading to steadily decreasing values of F. 

While the curve of F vs. S is not known exactly, for 
definiteness the portion of the curve between the 
origin and peak has been drawn to obey the following 
parabolic law: 


f = 2s — s? (2) 
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where, for convenience, the fractional force f and frac- 
tional slip s are referred to the maximum force and 
its associated value of slip respectively, thus: 


f = F/Fm = (3) 


QUALITATIVE DESCRIPTION 


As a preliminary, we shall examine the stability of 
the braking action on both the rising and falling parts 
of the curve of Fig. 2 with the spring assumed perfectly 
stiff, as clarified by Dr. J. V. Hughes and others, on 
the temporary assumption that the wheel inertia 
torque is negligible. On the rising part, if the angular 
velocity should for any reason fall slightly below the 
value corresponding to the torque 7 and its associated 
ground force F, then S will increase and cause F to 
increase, the latter tending to restore w to its equilib- 
rium value; thus, the action is stable. On the falling 
part, an incidental decrease in w and consequent in- 
crease in S leads to a decrease in F, thus decreasing w 
further, so that the action is unstable. These argu- 
ments are readily shown to apply also in the reverse 
sense. 

In practice, the falling part of the curve is normally 
entered after a stay on the rising part when, for any 
reason, the slip becomes greater than that required to 
reach the peak. With spring deflection and wheel in- 
ertia neglected, this can come about only through 
application of a braking torque higher than that corre- 
sponding to the peak of the curve. The excess of brak- 
ing torque over that required to equilibrate the ground 
force acts to decelerate the wheel, the deceleration in- 
creasing steadily as the slip rises to unity. When the 
wheel angular velocity finally reaches zero, then sliding 
of the brake surfaces stops and the torque falls instan- 
taneously to that required to equilibrate the ground 
force, so that rotational overshoot does not occur and 
the wheel remains locked. 

With the spring deflection taken into account, the 
action on both the rising and falling parts is much 
more complicated, but it is almost certain that violent 
disturbances, possibly including sustained oscillation, 
will occur if the falling part of the curve is reached. 
Even if the falling part of the curve is not entered, se- 
vere disturbances leading to gear walking may possibly 
occur if the system is uwmstable on the rising part at 
certain aircraft speeds. 

To explore these matters we first rewrite Eq. (1) so 
as to take ¥ into account (the dot denotes differentiation 
with respect to time ¢), thus: 


[Row/(V x) | (4) 


A typical maximum value of may be obtained from 
the experimental observation that a typical gear walk- 
ing frequency is of the order of 6 cps, with a double 
amplitude of the order of 6 in. This leads to a maxi- 
mum instantaneous value of x of about 7 m.p.h. Pre- 
sumably, the resulting disturbance in slip is unimpor- 
tant when V is many times this value but may be sig- 
nificant at low speed in producing instability or in raising 
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S beyond S,,. For example, if S, = 0.2 (a typical 
value) and the equilibrium value of S is 0.1, then for 
VY = 100 m.p.h. a value = —7 would only increase 
S to 0.16 (with w assumed constant), whereas for V = 
20 the slip would increase well beyond the value 0.2, 
presumably causing sudden locking of the wheel with 
resultant violent disturbances. 

Positive values of * are of no interest here because 
they tend to decrease S. If # > V, which can occur 
at very low values of I’, the wheel will overtake the 
ground, and S and F will be negative. 

The mass moment of inertia / of the wheel about its 
center cannot be neglected in the present discussion. If 
its value were zero, then F would always be fixed by 7 
so that S would be fixed, and the value of w would in- 
stantaneously be adjusted to suit x, according to Eq. 
(4). It is only as a result of the inertia torque /w# that 
F and, consequently, S can vary from their equilibrium 
values. 


MATHEMATICAL ANALYSIS 
Force and moment equilibrium in Fig. 1 require that 
kx + = F = (T+ Ia)/R (5) 


where R is the height of the wheel center above ground, 
and may be somewhat different from Ry in Eq. (1). It 
will be convenient, both for the purpose of simplifying 
the mathematics and aiding the intuition, to refer the 
various physical parameters to certain fixed values 
with the aid of new nondimensional or fractional 
parameters. Thus, we let x,, be the spring deflection 
when acted on by a static force F,,, and 2 the natural 
angular velocity (27 times the natural frequency) of 
the spring-mass system, as follows: 


Xm = F,,/k 
2 (6) 
so that, on introducing the fractional displacement 
= Xx/Xm (7) 
we have 
u + (u/2?) = f = (T/RF,,) + (Mr’Ro/ F,,) 
where = M(rR)*? (8) 


(r?AS,,/2(1 — s)]u’’ + u"[1 + 71 — sS,)] + u’f[PAS,, 
u = fo[l — (A 


Typical numerical values are thought to be as fol- 
lows: S, = 0.2; r? = 1/3, corresponding to a wheel 
mass of about 0.6.17 concentrated at a radius of gyration 
of about 0.75R; and fy a substantial fraction of unity— 
say, 0.75 (values of unity or over cause wheel locking 
under steady conditions). The value of s and its 
associated quantity f normally vary from zero in free 
rolling to a substantial fraction of unity. 

To estimate other numerical values, we make the 
rough assumption that oscillations in x are simple har- 


in which rR is the effective wheel radius at which the 
fictitious mass 7 would be concentrated to simulate 
the actual mass moment of inertia /. If we define the 
further dimensionless quantities 


time angle = 
fractional aircraft speed A = J Qy,, (9) 


fractional wheel speed w = Rw I" \ 
then we have, for the rising part of the curve of Fig. 2, 
u+u” =f = (T/RF,) + rAw’ = 2s s* (10) 
SSm = 1 — [w/(1 — u’/A)] (11) 


where use has been made of Eqs. (2) and (4); the 
former being only a rough approximation and the latter 
being approximated by replacing Ry with R; and 
primes denote differentiation with respect to ¢. 

We must now allow for an essential factor in the 
analysis—namely, the variation of the torque 7 with 
angular velocity w. Conventional data* on automobile 
brakes indicate that 7 decreases gradually with increas- 
ing relative velocity of the rubbing parts; over a rea- 
sonably small range of w, therefore, such as would be 
covered if the slip never exceeded the value S,, in Fig. 
2, we shall adopt the following approximation: 


T = — w/wy) (12) 


where 7) is the fictitious value that would be reached 
at w = Oif the linear relationship were extended to this 
value, and, similarly, w, is the fictitious value of w at 
which 7 would disappear. We shall use this equation 
with the reservation that for our present purpose it is 
not adequately supported by test data, since in con- 
ventional frictional experiments no particular attention 
is paid to the time factor, whereas in the present investi- 
gation the behavior for fast oscillations of w (of the 
order of 6 per sec.) is of interest. 
By Eq. (11) and the third of Eq. (9), we have 


T = T){1 — (V/Rey) [1 — (u’/A)] (1 — sS,,)} 


Rw; may be interpreted as an equivalent aircraft 
speed V’,, represented by a corresponding relatively 
high value of A—say, A,—so that the term in the first 
parentheses may be written as A/A,. Writing 7)/R 
as F, and F,/F,, as fo, we then have, by Eqs. (10) and 
(19): 


2(1 — s)] — [fo(1 — sS,)/Ay]} + 
A,) (1 — sSn)] + [r?Sp/2(1 — s)]u’(u’ + (13) 


monic, with any given amplitude x») and the natural 
angular velocity 2 of the spring-mass system, thus: 


x = X sin Qt (14) 


| 


Then 
Xmazr = X0, Xmazr = Xmazr = 


so that, by Eq. (7) and the first part of Eq. (9), u and 
its derivatives with respect to ¢ are of the order of 


frac- 
e and 
| 


unity when x) = x», that is, under very severe vibration. 
(This convenient result is a consequence of the choice 
of dimensionless parameters.) To find significant 
numerical values for A, we assume that {2 corresponds 
to a typical observed gear walking frequency of 6 eps, 
and that x,,, the landing gear «-deflection under a static 
load F,,, has the typical value of 8 in., so that the value 
A = 1 corresponds to an aircraft speed of about 17 
m.p.h. 

While Eq. (13) is complicated, the stability for small 
disturbances about a steady value of u can be investi- 
gated readily. As the steady value (with all deriv- 
atives of u equal to 0) take 


+ [201 — 1)/r?AS,, + 


where terms of lower order of magnitude have been 
dropped. 

By the numerical considerations discussed previously, 
plus the requirement that Ay be very much larger 
than A, the final term in the braces multiplying 6” 
can be ignored in comparison with the remainder of the 
terms in the braces, and, similarly, for the braces 
multiplying 6. Also, the term S,,s; can be replaced with 
little error—in the two places where it appears—by the 
typical value 0.2 X 0.5 = 0.1. Then the above equa- 
tion reduces to: 


+ [2(1 — s,)/r?AS,, ]6"(1 + 0.97?) + 
— [1L.8(1 — si)fo/r?AS,A + 
— s1)/r?AS,,]6 = 0 (18) 


By conventional stability theory, the system will be 
stable against small disturbances if none of the coeffi- 
cients of 6 or its derivatives falls below zero, or if the 
product of the second and third coefficients is not less 
than the fourth. If A, is the critical value of A, namely 
that at which the stability becomes marginal, it can be 
seen by inspection that the second of these two criteria 
dominates since it represents a higher value of A,, below 
which all values of A represent instability. Applica- 
tion of the second criterion yields: 


A, = 2(1 — s1)fo (1 + 0.9 r?)/r4S,,A 7 (19) 


where, in view of the fact that A/A,in Eq. (15) is very 
small when compared to unity, fo is closely equal to 
251 — 51”. 

Let us assume the value A; = 30, corresponding 
(with previous numerical values) to a fictitious aircraft 
speed of about 510 m.p.h. at which the braking torque 
would theoretically disappear. (This value of A, is, 
as referred to the relatively low aircraft speeds at which 
instability is observed, extremely conservative in the 
light of the previously cited data of Gemant.*) Then 
instability is indicated when A, = 1.5—that is, when 
the aircraft speed is about 26 m.p.h. or less. Because 
of structural damping, neglected in the above analysis, 
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where Eqs. (13) and (10) have been used, and 


Ss Si t 


(16) 


/ 


where, by the previous numerical discussion, 6 and its 
derivatives, as well as e, are very small compared to 
unity. Using Eq. (16) in conjunction with the first 
and last of expressions in Eq. (10), we find: 


e = [(6 + 6")/2(1 — (17) 


Substituting Eqs. (15)-(17) in Eq. (13), and taking 
general note of numerical values, including the fact 
that values of A much smaller than unity represent air- 
craft speeds too small to be of practical interest, we find: 


SnS1) Snfo 2A (1 si) |} + 1 [2(1 fo 7] 
(1 — Snsi)} + [2 — 6{1 — [A Spfo/24,1 — = 0 


instability may actually be delayed until the aircraft 
speed has fallen to a somewhat lower value. 

For easy interpretation of the above equation, we 
shall change certain of its fractional parameters to 
actual physical quantities. First we rewrite Eq. (8) 
thus: 

M,R,? = I =(M, + M,) (rR)? 
where J, is the wheel mass, 1/, is the remainder of \f 
and consists of part of the mass of the landing gear and 
its connections, and R, is the radius of gyration of the 
wheel. Substituting the resulting value of rin Eq. (19), 
and writing V. and V’,;in place of A, and A, with the 
help of the second part of Eq. (9) and also Eq. (6), we 
get: 
Ve = 2(1 — Si)fotmPm[M, + My (1 + 0.9 7,7)] + 
(My 1,7)? (20) 
where 7, is the ratio R,/R. To eliminate V,, we note 
the quantity 07'/Ow, which has been assumed constant 
and negative. Since V, is a small fraction of V,, and 
since the slip is small so that Rw is roughly equal to V, 
we may write the following approximation : 


To/Vs = —[0T/0(Rw) | 
where 7) is the braking torque corresponding to the 
static fractional ground force fy). Also, we define the 
following fictitious deceleration : 

a = Fy/{(Mur,?)?/(M, + Mw (1 + 0.9 r,2)]} (21) 
so that Eq. (20) finally becomes 
2(1 — 51) fodxm | - | 


Ve 
Sin 


(22) 


where the partial derivative has been written in a form 
in which it is relatively insensitive to the greatly varying 
values of fo. 

There is little prospect of obtaining significant relief 
from gear walking through the modification of any of 
the parameters in this equation, except for the partial 
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derivative, which may possibly be reduced or, better 
yet, made positive by suitable braking materials. As 
to the term (1 — s;)f) which is fixed by the brake pedal 
force exerted by the pilot, this must be left free to vary 
over a wide range to suit different braking situations. 
The value of a is—if we assume that 7, cannot be 
greatly varied and that J/,, M,, and the overall mass 
of the aircraft are proportionately reasonably fixed by 
general design considerations—indicative of the maxi- 
mum aircraft deceleration during the landing run, and 
so cannot be varied much. It is, however, of theo- 
retical interest to note by Eq. (21) that a decrease of 
M, and an increase of M/, would be favorable. The 
value of x, = F,,/k has presumably already been made 
as small as possible by maximizing the stiffness k. 
Finally, S,, is dependent on the properties of the tire 
and the landing field, which, presumably, cannot be 
intentionally varied in such a way as to increase S,, 
significantly. 


Sinusoidal Analysis of Zero-Stability Condition 


A detailed physical understanding of the conditions 
at which the small disturbance stability just vanishes 
can be gained by assuming that a small sinusoidal dis- 
turbance is imposed on the system by any external 
means, and then noting the requirements for the dis- 
turbance to be sustained indefinitely after the external 
agency is withdrawn. Hence, assume 


Au U sing ¢ \ = 
(23) 
Aw = W sin (qg — a)f 
where LU’ and W are small amplitudes, g is the ratio of 
the angular velocity of the disturbance to the natural 
angular velocity 2 of the spring-mass system, and a is 
an angle of phase lag. U and g are prescribed, while 
W and a are unknown. However, the value of g at 
which the oscillations are self-sustaining will probably 
turn out to be unique, and, by intuitive considerations, 
will be slightly less than unity because of the moment 
of inertia of the wheel. 

Introducing Eq. (23) in Eq. (10), making use of 
Eq. (12) and the second and third part of Eq. (9), and 
subtracting out terms relating to the steady state of 
motion, we find: 


— q*) sin gg = Af = —fo(A./A)W X 
sin (gg — a) + r°A.gW cos (qg — a) = 
2(1 — s;) As (24) 


where Af and As are small disturbances resulting from 
the applied disturbance in uw, and terms of higher order 
in As have been dropped. The subscript c has been 
written for A since we are limiting the present investi- 
gation to the zero-stability condition. 

By similar considerations we may eliminate As with 
the aid of Eq. (11), obtaining: 


Af = [2(1 — s:)/S,] [—W sin (qe — a) — 
(1 — (qU/A,) cos (25) 


Fig. 3 shows a conventional time-vector diagram 


Ac 


- W sin(q@-a) sing@ 


PA. qw cos(q@-a) 


Af=U(I-q*) sin q@ 


2 (I-s,)(I-s, S,) 
SA. qU cos 


a(i-s,) A, 
Ay W sin(q?-d) 


Aw 


A, q W cos 


sin 


Fic. 3. Time-vector diagram. 


which satisfies Eqs. (24) and (25). The triangles have 
been drawn roughly in proportion, in accordance with 
numerical quantities discussed previously. Making 
certain obvious approximations, and noting especially 
that A, < Ay, we find: 


fo(Ac/Ay) [201 — 81)/Sm] = 
— = rAqw 


(26) 
W= fe! A,|qU 


These lead to 
g = 1/[1 + r(1 — s,S,)] 


so that gq is slightly less than unity, as expected. Solv- 
ing for A, and approximating s,5,, by 0.1, as before, we 
find Eq. (19) is confirmed. 

It is of interest to examine the instantaneous powers 
P, and P, absorbed at the brake and ground, respec- 
tively. Using previous definitions freely, we may 
write: 


P, = Tw = F,,R Aw] x 
(V/R) (w: + Aw) 
P, = —F(x + Rw) = —Fy(fi + Af) X 


[u'(V/A.) + + Aw)] 


Defining p, and p, as the respective ratios of power, 
other than that corresponding to the steady motion, to 
F.,, V, we have: 


= fidw A pw Aw = ‘A (Aw)? 
pb, = —filu’/A) — fidw — wAf — (u'/A Af — Afdw 


All terms consisting of the product of a constant by 
a sinusoid vanish when integrated over a cycle of 
vibration. The next to last term in the second equation 
also vanishes when thus integrated because in Fig. 3 u’ 
leads Au by 90° and thus is in quadrature with Af. As 
to the last term in the first equation, this is always 
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negative; if we write the last term in the second equa- 
tion as follows with the help of the second and third 
expressions of Eq. (24): 


[fo(A./A — r?A.(dAw/dg) |Aw 


then the last terms in the two respective equations are 
seen to have equal and opposite integrals over the 
cycle, so that energy is transferred from brake to ground. 

That the brake has a negative increment of energy is 
not surprising in view of the fact that the assumed 
negative slope of its curve of torque vs. angular ve- 
locity gives it the characteristic of a negative dashpot. 
Conversely, the surface between tire and ground must 
have the characteristic of a positive dashpot, since a 
positive increment of slip, which produces a positive 
increment of ground force, is associated with an in- 
crease in relative velocity of the wheel in the direction 
opposed to ground force. 


Moderate-Disturbance Stability 


In the above stability analyses, the curve of braking 
torque T vs. angular velocity w has been assumed free 


Am + — q’) singe + U2(1 — 4g’) sin — 8) = 


of curvature, and the curvature of the curve of ground 
F vs. slip S has been ignored, so that the results apply 
only for disturbances which are sufficiently small so 
that the second-derivative terms in the Taylor ex- 
pansion for the curves are negligible. We next con- 
sider somewhat larger disturbances—namely, of such 
a magnitude that the second-derivative terms are no 
longer negligible although still relatively small when 
compared to the first-derivative terms. With this 
assumption, it will be shown that a disturbance will be 
self-sustaining if Eq. (23) is replaced by the following: 


Au = Am + Usin gg + U2 sin — 8) 
Aw = Aw, + W sin (q¢g — a) + > (27) 
where the additional unknown parameters are the 
steady increments Aw, and Aw,, the second-harmonic 
amplitudes l’, and W2, and the phase lag angles 8 and 
y. Applying these to Eq. (10), as before, and taking 
due note of our assumptions as to the magnitude of 
various terms and of the fact that A, < Ay, we get 
in place of Eqs. (24) and (25): 


Af = —fo(A-/Ay)W sin (gg — a) + 


fo(b/4)W? [1 — cos — a)|] + [W cos (gg — a) + 2W2 cos ¢ — = [—2(1 — 51)/Sn] X 
{ Aw + W sin (qg — a) + We sin (2g — vy) + (wig/A.) [U cos gg + 2U2 cos (2g¢ — B)| + 
(qUW/2A,) [sin — a) — sin a] + (wig?U?/2A.2) (1 + cos + [(—1)/Sn?] X 

{(W? 2) [L — cos 2(gg — a)| + (1/2) (wigU/A,)? (1 + cos + (wigUW/A.) [sin (2g¢ — a) — sin a]} (28) 


where b is the second derivative of 7/7, with respect 
to w, and the second derivative of /, with respect to s, 
namely —1, is displayed in the parentheses before the 
last pair of braces. 

These equations must be separately satisfied by the 
constant, first-harmonic and second-harmonic terms. 
The first-harmonic terms are seen to be identical to 
those contained in Eqs. (24) and (25), if we make use 
of the fact that w; = 1 — s,S,,.. To solve for the con- 
stant and second-harmonic unknowns it will be con- 
venient to adopt the following specific numerical 
values, which are consistent with numerical values dis- 
cussed previously: 

fo A; U 
0.75 0.5 0.2 30 1/3 0.2 
and where a substantial value has been chosen for L’. 
Eq. (26) and Fig. 3 then yield: 

q A, W a 
0.877 1.46 0.108 95° 
where the value of a justifies the narrowness of the tri- 
angles in Fig. 3 and the resulting approximations 
leading to Eq. (26). To apply these values to Eq. 
(28), we assume that 6 = 1, which can be shown to be 


far greater than can be expected in practice, with the 
following results: 


Au, Aw, U2 B 
0.0022 —0.00044 0.0007 0.0004 108° 173° 


where the magnitudes of the various terms justify the 
basic numerical assumption stated at the beginning of 
this section. 

In summary of this portion of the work—for moder- 
ate disturbances where the terms arising from the curv- 
ature of the F vs. S and T vs. w curves are appreciable 
but still relatively small, the conditions for stability 
are the same as for very small disturbances. For 
sinusoidal disturbances in particular, relatively small 
constant and second-harmonic increments are super- 
imposed on the first-harmonic terms. 


Gear Walking Due To Braking Torque Overshoot 


Even though the slope of the curve of braking torque 
vs. angular velocity is not negative enough to result in 
instability, it is possible for severe vibration of the 
landing gear to result from sudden application of the 
brake. To explore this, consider the particular case 
in which the aforesaid slope is zero, so that, with the 
typical numerical values mentioned previously, Eqs. 
(13) and (18) reduce to the following: 


[4/3001 — s)] + 1.30” + [A/30(1 — s) Ju’ + 
= fo +[u'(u’ + u’’’)/30(1 — s)] (29) 


+ [30(1 — + 8’ + 
(301 — s:)/A] =0 (30) 


By previous criteria, the latter equation shows the 
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small-disturbance stability to be always positive. 
However, at low values of A, with consequent high 
values of the coefficients in the second and fourth 
terms, we may expect the damping represented by 6’ 
to be relatively unimportant, with the result that the 
overshoot in “, and, consequently, in f, resulting from 
sudden brake application, may be large enough to cause 
f to exceed unity instantaneously, with consequent 
locking of the wheel, as explained in the chapter— 
Qualitative Description. To examine this for definite 
values of A, it was assumed that the wheel was rolling 
freely and the value f/f) was then applied suddenly. 
The step by step solutions of Eq. (29) for the resulting 
behavior were found for A equal to 15 and 3, correspond- 
ing (with previous numerical values) to aircraft speeds 
of about 250 and 50 m.p.h., respectively, and are 
shown in Fig. 4 as plots of f vs. ¢. 

As expected, the overshoot of f beyond its equilib- 
rium value is seen to be negligible in the high-speed 
case, but quite appreciable in the lower-speed case, 
though not enough to reach the locking value of unity. 
Undoubtedly, this would have been reached in the 
latter case if fy had been chosen somewhat higher. If 
this happens in practice, then the brake pedal is likely 
to be released completely to allow the wheel to unlock 
and to roll freely, followed by sudden re-application 
of a force corresponding to fy, thus leading to the same 
overshoot as before; if this cycle is repeated several 
times, very severe oscillations result. 


REMEDIAL MEASURES 


As shown previously, very little can be done in a 
critically designed aircraft to modify the parameters, 
other than the slope of braking torque vs. angular 
velocity, so as to secure significant relief from gear 
walking. If research fails to turn up satisfactory 
braking materials having negligible or positive slope, 
then recourse may be taken to a damping device to 
supply either x-damping or w-damping. 

For x-damping, a dashpot may be connected between 
the landing gear and relatively fixed portions of the 
aircraft, or possibly the landing gear may be furnished 
with loose parts (conceivably in the form of pellets in 
hollow spaces of the gear) which collide with or rub 
against each other during oscillations, and thus absorb 
energy. For w-damping part of the braking may be 
done by a fluid brake in which the torque increases 
with increasing angular velocity, or else, some form of 
rotary inertia damper may be fitted to the wheel. The 
fluid brake must be capable of withstanding the very 
high angular velocity at the beginning of the landing 
run, and the inertia damper must not lead to excessive 
wear of the tire at the instant of touchdown. 

With an appreciably positive slope, gear walking 
will not take place unless excessive brake-pedal force 
is repeatedly applied. If the slope is negligible, then 


the system has been shown to be stable for small de- 
partures from equilibrium conditions, but there is ob- 
tained a strong overshoot in f for low values of relative 


A=15 


| 2 3 


4 5 
l l 


Fic. 4. Fractional force f vs. time angle ¢. 


aircraft speed A when the initial conditions are those 
resulting from release and subsequent re-application 
of braking torque after locking of the wheels has 
inadvertently taken place. Hence, a straightforward 
remedial measure is to decrease the braking torque and, 
therefore, fy, as the aircraft slows down, so that / never 
exceeds unity even with a large overshoot. Since the 
curve of Fig. 2 changes for different conditions of normal 
wheel force and landing field surface, these variables 
must be taken into account to obtain satisfactory values 
of fo. At lower speeds, it is also important, in order to 
prevent appreciable overshoot, to apply the brake cau- 
tiously after it has been fully released, as, for example, 
after inadvertent wheel locking. Systematic methods 
for accomplishing this will not be described in detail 
here. The damping devices discussed above will, in 
general, also serve to reduce overshoot. 


SUGGESTED EXPERIMENTS 


To apply the above theory, it is essential that the 
slope of the curve of braking torque vs. wheel angular 
velocity, for rapid variations of angular velocity (changes 
of the order of 10 r.p.s. per sec.) be determined for 
braking systems of actual interest in connection with 
gear walking. It may or may not be necessary to in- 
clude oscillations in the variation in order to get indica- 
tive data. 

If the slope of the curve is found to be negative, then 
the damping devices discussed above might be investi- 
gated experimentally, either on a full-scale airplane or 
on a model corresponding to Fig. 1. If, on the other 
hand, the slope is zero or positive, then the experi- 
ments should aim at determining whether gear walking 
can be avoided merely by avoiding large initial dis- 
turbances in braking and by keeping the braking 
torque well below the locking value. 

Additional theoretical and experimental data on the 
shape of the rising part of the curve of ground force 
vs. slip would be useful in checking the suitability of the 
parabolic approximation used herein. 

(Continued on page 320) 
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Thermal Stresses in Rings 


Marvin Forray 

Design Specialist, Republic Aviation Corporation, 
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November 24, 1958 


HE PROBLEM Of a traction-free solid plate subject to a general 
two-dimensional temperature distribution was considered in 
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)) is 


(2) 


(3) 


op = gor) + [gn(r) cos + h,(r) sin (4) 


n=1 


where the functions g,, h, are obtained from 
yan S S dr dr (n = 0,1, 
Sr S " dr dr i233...) f 


The stress components are obtained from @ by the formulas 


(5) 


Orr = (1/r) (0¢/dr) + (1/r?) (0%6/06?) 
= —(0/dr) [(1/r) (0¢6/08)] > (6) 
o99 = (0°6/Or?) 


CrRCULAR RING 


Consider a circular ring of inner radius }, outer radius a, sub- 
ject to a temperature distribution 


T = (r/a)? (1 — cos 6) + (7) 
where 79, 7; are constants and the polar axis (@ = 0) is assumed 
to be vertically upward. Thus Py) = —Ea[7T; + (r/a)?*]; 
P, = EaT((r?/a?); = 0,n >1; Qn =0,n =1,2,... Substi- 


tution of the quantities P,, Q, into Eqs. (5) and (4) yields 
op = (—Ea/4) [Tir? + (Tort/4a?)] + (EaTor*/15a?) cos (8) 


In order to ensure single valuedness of the stress and displace 
ment components? for traction-free boundaries 


= dy = ay" =a, = b,’ = = d,’ = 0 (9) 


The satisfaction of the traction-free boundary conditions requires 
that ¢- assume the reduced form 


db. = dy log r + bor? + (bir? + a,’r~) cos 6 (10) 
Substitution of Eqs. (10) and (8) into Eq. (6) yields 
Orr = (ao/r?) + [2b9 — (EaT;/2)] — (EaTor?/4a?) + 
{(EaTor?/5a?) + 2bir — (2a;’/r*)] cos 6 
[((EaTor?/5a2) + 2hir — (2a,’/r*)] sin 6 (11) 


—(ao/re) + [2b9 — (EaT,/2)] — (8EaTor?/4a*) + 
+ + (2a;'/r*)] cos 


a 
> 
| 
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The coefficients ao, bo, b:, a; are determined from the boundary 


conditions = = Qonr =b: 


ado = —(EaT)/4)b? 
by = (Ea/4) {7, + (To/2) [1 + (b?/a?)]} (12) 
a,’ = —(EaT)/10) [a*b4/(a + b) (a? + 
by = —(EaT)/10a) {1 + [b*/a(a + b) (a? + b?)]} 
so that the stress components in nondimensional form are 
1 (- b? 414 
EaT»o 4 r? a? 
5 la? a a(a + b) (a? + b?) 
a*b4 
(a ra b) (a? + b2)r3f 
EaTy 5 (a? a a(a + bd) (a? + 5b?) 
(13) 
a*h4 


(a b) (a? + 
a? 
cos far? 3r b4 
5 Va a a(a + b) (a? + 6?) 


a*h4 
(a + b) (a? + 


which is readily verified to be the correct solution. 
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On the Representation of Nonlinear Creep by 
a Linear Viscoelastic Model* 


Harry H. Hilton 
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INEAR VISCOELASTIC STRESS-STRAIN LAWS have found much 
favor in high-temperature structural analysis because mathe- 
matically they are relatively simple to handle and because they 
lend themselves to analogies with equivalent elastic problems.'~* 
However, in stability problems, such as creep buckling and creep 
divergence, linear viscoelastic studies have not produced the 
required critical times,‘ > while the use of nonlinear viscoelas- 
ticity has, with a few exceptions,® ® lead to finite times. 

The experiments of Allen, Brull, and Wilkie* have shown that 
linear viscoelastic models such as the standard linear solid and 
Maxwell-Kelvin solid fit well the experimentally determined 
creep curves under constant load and temperature for 2017 and 
2024 aluminum alloys. Baer’ has independently derived non- 
linear expressions for creep strains from experimental data on 
QQ-M-44 magnesium alloys. 

It will be shown here that these two results, instead of being 
contradictory, actually are identical. That is, it will be proved 
that under constant temperature and load the Maxwell-Kelvin 
model corresponds to Baer’s empirical creep formula and, conse- 

* The analysis described in this paper was performed under the sponsor- 
ship of the Systems Development Laboratories, Hughes Aircraft Company, 
and the author is grateful to them for permission to publish this article. 
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quently, under these conditions the behavior of these metals is 
the same. 

The Maxwell-Kelvin body may be represented by a mechanical 
model consisting of two linear springs and two linear dashpots, 
arranged as shown in Fig. 1. The stress-strain relation in this 
model is for uniaxial stress 
(0/dt)[3e — (@/K) — 3aT] + (0/dt) 7(0/dt) X 

[Be — (@/K) — 8aT]} = (0/dt)[72(0/dt) + 
(0/ot) + o} + (a/m) (1) 


where rr = m[(1/Gi) + (1/G2)] (2) 
n= no/G» (3) 


and with K, a, and 7, respectively, the bulk modulus, coefficient 
of linear thermal expansion, and the temperature. The initial 


conditions are 
3(e — aT) = of[(1/K) + (1/G;)]  (t = 0) (4) 


3(0/dt) (e — aT) = (0/dt) + + 
a((1/m) + = 0) (5) 


For the case under consideration, where both temperature and 
stress are time-independent, Eq. (1) reduces to 


+ = o/3m (6) 
The solution of Eq. (6) subject to the conditions (4) and (5) is 
= +(1/K)] + (t/m) + — (7) 
Comparing this result to Baer’s empirically determined equa- 
tion 
= Co™ + + C1 — (8) 
shows that one need only to select the model parameters such that 


= — (Y/R) 


I/m = 3C.0"~! 
1/72 = C; 


to complete the correspondence between the nonlinear creep law 
(8) and the Maxwell-Kelvin model. 

It should, of course, be recognized that the foregoing applies 
only to time-independent temperatures and stresses. If the tem- 
perature and/or stresses are time functions, then Eq. (1) rather 
than Eq. (6) describes the linear model and, furthermore, the 
validity of Eq. (8) under variable conditions has not been demon- 
strated. 

It is also of interest to note that for a viscoelastic body, such 
as Maxwell-Kelvin model, it is not possible to predict variable 
stress creep from constant stress and temperature experiments. 
First of all, Eqs. (1) and (6), describing these two phenomena, 
are radically different. Secondly, the number of material con- 
stants describing the variable conditions are five (Gi, K, 71, 72, 1), 
while only four are needed for constant stress conditions. 
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Further Note on Isentropic Compressible Flow 
of Air With Variable Specific Heats 


K. E. Tempelmeyer 
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December 1, 1958 


 yrewogons FLOW EQUATIONS which included real-gas effects 
were derived for air at temperatures up to 4,500°R. in refer- 
ence 1. The effects of specific heat variation with temperature 
were considered by subdividing the c,-7 curve into three 
ranges and empirically fitting an equation into each range. 
These empirical relationships were then employed in the integra- 
tion of Euler’s and the general energy equation. Several sets of 
equations result; these depend upon the ranges in which the 
static and stagnation temperatures lie. Because of space limita- 
tions, only the derivation was presented in the original note. 

The usefulness of these equations in the range between perfect- 
gas compressible flow relationships and real-gas hypersonic flow 
charts has been pointed out in reference 2. The complete set of 
these equations for isentropic flow is tabulated below in order to 
facilitate their use for calculations in this intermediate range. 
The notation employed is identical to that of reference 1. All 
empirical constants from the specific heat equations have been 
inserted and simplified. The tabulation below is intended as 
an appendix for reference 1. 


(1) 7,and in Range 1 
(1) cp = 0.2393 
(2) = ¢p/[ep — (R/J)] 
(3) M? = [(2/) (0.2393) — T)]/yRT 
(4) p/pe = (T/T1)** 
(5) = (T/T,)?™ 


(II) 7,and T in Range 2 
(1) cp = 0.2318 + 0.1040 X 10-4T + 0.7166 XK 10-®T? 
(2) = ¢p/[ep — (R/J)] 
(3) M? = 2J[0.2318(T, — T) + 0.0520 K 10-4(7;? — T?) + 
0.2388 X 10-8 — T?)]/yRT 
(4) = (T/T,)**™ exp [(1.516 X 10-4) (T — T,) + 
(5.224 10-5) (T? — T,*)] 
(5) p/p: = (T/T,)?*™ exp[(1.516 XK (T — T1) + 
(5.224 10-8) (T? — T,?)] 
(III) 7,and 7 in Range 3 
(1) cp = 0.2214 + 0.3521 X 10-*T — 0.3776 10-87? 
(2) vy = ¢p/[ep — (R/J)] 
(3) M2? = 2/[0.2214(7, — T) + 0.1760 X 10-“ 72 — T*) — 
0.1258 — T%)]/yRT 
(4) p/pe = exp[(5.134 X 10-4) (T — 
(2.753 X 10-8) (7? — T,?)| 


(5) p/pr = (T/T ,)2229 exp [(5.134 x 10 4) (T = T,) 
2.753 X 10-8) (T? — T,*)] 


(IV) T,in Range 2 and T in Range 1 
(1) cp = 0.2393 
(2) = ¢p/[ep — (R/J)] 


(3) M? = (2/) [(0.2214)T, + (0.0520 X 10~4)T,? + 
(0.2389 X 10-8)T,3 — (0.2393)T]/yRT 


(4) = (T/400)* (400/T,)* exp[0.06901 — 
(1.516 X 1074)T, — (5.224 + 
(5) p/p: = (T/400)249" (400/T,)2™ X exp[0.06901 — 


(1.516 X 10~4)T, — (5.224 X 10-°)T,?] 
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(V) T;1n Range 3 and JT in Range 2 
(1) cp = 0.2318 + 0.1040 K 10~4T + 0.7166 XK 10-87? 
(2) = — (R/J)] 


(3) M? = (2/7) [(0.22147, + 01760 10-*T/? — 0.1258 X 
10-873) — (0.2318T + 0.0520 X 10-472 + 


0.2388 10-8T3)|/yRT 


(4) p/p: in (T/1700)?-389 > 4 
exp[(2.753 XK — (5.134 X 10~*)T, + 0.7932] X 
exp[(1.516 X 107*)T + (5.224 X 10-8)T? — 0.4087] 
(5) = (T/1700)?-™ (1700/T,)2-229 X 
exp[(2.753 X — (5.1384 X 1074)7T; + 0.7932] X 
exp[(1.516 10°*)T + (5.224 1078)T? — 0.4087] 
(VI) 7, in Range 3 and 7 in Range 1 
(1) cp = 0.2393 
(2) = ¢p/[ep — (R/J)] 


(3) M? = (27) [(0.22147, + 0.1760 X — 


0.1258 10°-8T,3) — 0.23937] /yRT 


(4) = (7 /400)*-49! XK (0.00535) X 
exp[(2.753 — (5.134 & 1074)T, + 0.7932] 


(5) p/pe = K (T/400)?-#! (0.0228) X 
exp[(2.753 10-8)T? — (5.134 X 10~4)T, + 0.7932] 
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On Flat Plates and Detached Shock Waves 


S. S. Sopezak 

Aerodynamics Department, Goodyear Aircraft Corporation, 
Akron, Ohio 

December 5, 1958 


ged UNDERGRADUATE THESIS work by the author at the 
University of Detroit leads to the possible modification of 
Serbin’s Eq. (5) of reference 2. This modification would make 
possible the prediction of the shock detachment distance for any 
flat-plate configuration. As stated in reference 2, the shock 
detachment distance for a disc can be determined by the following 


formula: 
A/R = 1.03/(K — 1)!/2 (disc) 


where A is the detachment distance, R is the radius of the disc, 
and K is the density ratio across a normal shock at the free- 
stream Mach Number. 

The modified formula is: 


A = 0.581 [A/(K — 1)]!/" (flat plate) 


where A is the area of the flat plate. This modified formula was 
verified by the test results obtained on flat plates of a circular, 
square, and rectangular configuration. Tiie tests were limited 
toa Mach Number of 2.82, and the aspect ratios of the flat plates 


were less than 1.5. 
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On the Flow of a Hydromagnetic Fluid Near an 
Oscillating Flat Plate 


R. S. Ong and J. A. Nicholls 

Department of Aeronautical Engineering, University of Michigan, 
Ann Arbor, Mich. 

December 8, 1958 


bw IMPULSIVE MOTION of an infinite flat plate in a viscous, 
incompressible magnetic fluid in the presence of an external 
magnetic field has been discussed by Rossow.'! In this note, 
his method is extended to cover the case of the flow near an in- 
finite flat wall which executes linear harmonic oscillations parallel 
to itself. The velocity profile will be found for the two cases: 
(1) the magnetic lines of force fixed relative to the fluid, (2) the 
magnetic lines of force fixed relative to the plate. 


MAGNETIC FIELD FIXED RELATIVE TO THE FLUID 
Let x denote the coordinate parallel to the direction of motion 
and y the coordinate perpendicular to the wall. At time ¢ < 0, 
the fluid, plate, and magnetic field are assumed to be stationary 
everywhere. The plate starts to move at time ¢ = 0. 
The flow may be approximately described by the following 
differential equation 

(Ou/dt) + u = v(02u/dy?) (1) 
where: « = x — component of fluid velocity; By = external 
magnetic field directed perpendicular to the plate; « = con- 
ductivity of the fluid; » = kinematic viscosity of the fluid. 
Bo, o, and vy are assumed to be constant. The boundary condi- 
tions are: u = ujcosntat y = 0,t > O;and u = finiteat y = 
t>0. 

The Laplace transform of the velocity u is defined as 


i= e~*t udt 
0 


Applying the Laplace transform to Eq. (1), gives 

v(d*i/dy?) — (s + m)a = (2) 
where m = aB,?/p. 

The solution of Eq. (2) is given by 

ti = C\(s) exp[—y{(s + m)/v}! 2) 4 CAs) expLy}(s + m)/v}! ‘2) 
Since u = finite at y = ©, C.(s) is chosen to be zero. The inte- 
gration constant C,(s) is found from the boundary condition at 
y = 0—i.e., on the plate—giving 


Ci(s) = uos/(s? + n?) 
Hence, 
ii(y, s) = {uos/(s? + n®)} exp [—y](s + m)/v}!/7] (3) 


Eq. (3) may now be inverted by means of the inversion theorem :? 


50 | —- -4 
| a vs 
40. 
B=O0, nt=O 
B O, nt=m7/2 
3! 3. B=10,nt=0 | 
a | 4 B=10, nt=7/2 
2 __| 5. nt=0 
LS. B= 10, nt=72, | 
| | | 
| 
-0O4 -O2 O 04 06 08 LO 
Fic. 1. Velocity profile with magnetic field fixed relative to the 
fluid 
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a_vs. 

4 B=0, nt=O 
2.B=0, nt=7/2 
3. B= 025, nt =O 
\ 4. B= 025, nt=7/2 


5. B= 1.0, nt=O 
6. B= 10, nt=7/2 


-02 o2 04 O06 O8 10 


u/Uo 
Fic. 2. Velocity profile with magnetic field fixed relative to the 
plate. 
y+iB 
u(t) = {1/27} lim ds (4) 
B- @ 


After a straightforward but lengthy computation, the solution 


muo 


ws, y) = (s + m) st +n? 


From the boundary conditions: C.(s) = 0, and C\(s) = s/m. 


Hence, 
(s + m) (s? + n*) 
s + m\'/?\ 
(s + m) (s? + n2) ) 
mM Uo 
u(y, t) m? n 4 n? 


Fig. 2 illustrates this motion for various strengths of the external 
magnetic field and for several instants of time. Again a = 
y(n/2v)'/2 and B= m/n. 

In both figures, the curves for mt = a and 37/2 are not in- 
cluded due to the complete symmetry (Fig. 1) and almost com- 
plete symmetry (Fig. 2) with respect to the ordinate. 
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Generalized-Newtonian Theory 
E. S. Love 


Head, Hypersonic Tunnels Branch, 
Langley Research Center, NASA, Langley Field, Va. 
January 6, 1959 


| Byes has shown that in application to blunt-nose bodies an 
improvement upon Newtonian theory is realized by use of 


MAY, 1959 


to Eq. (1) may be expressed by 


u(y, t) = uo exp { —y[(1/2r)((m? + n?)!/2 + m)]! 3} x 
cos {nt — y[(1/2v)((m? + n?)!/2 — m)}!/?} (5) 


It is to be noted that if m = O0—i.e., with no external magnetic 
field—the above expression reduces to the well-known solution 
for the ordinary viscous fluid flow problem.* Fig. 1 represents 
qualitatively this motion for several instants of time and for var- 
ious strengths of the external magnetic field where a = y(n/2v)'”? 
and 8B = m/n. 


MAGNETIC FIELD FIXED RELATIVE TO THE PLATE 


In this case, the magnetic field is moving and the fluid is 
initially at rest, so the relative motion must be accounted for. 
Hence, the relation analogous to Eq. (1) is 


(Ou/Ot) + m(u — uo cos nt) = v(d2u/dy?) (6) 


The boundary conditions are: u = um cos nt at y = 0,t > 0; 
and u = finite at y = ~,¢ > 0. Application of the Laplace 
transform to Eq. (6) gives 


/dy? = (m + — muo{s/(s? + n?)} (7) 


The solution of Eq. (7) is given by 


+. + m) | Gus) exp -y ( ) + exp {y( ) f (8) 


The first term on the right-hand side may be inverted by ele- 
mentary methods to give 


( 


(s + m) (s? + n?)f 
m? + 
The second term is inverted by means of the inversion formula, 
Eq. (4), and the final solution to Eq. (6) is expressed by 


E cos nt + n sin nt — m exp (—mt | 


‘[m cos nt + n sin nt] + ———— exp4-— y (m2 + n2)/2 4+ om t x 


1 1/2 1/2 
[» cos -y E (om n2)!/2 m)| — msin -y E G + n2)l/2 — n) | (10) 


the modified form 
C5/Cy = sin? 6 (1) 


max 


that is derived from the normal 


together with the value of C), 
shock relations. (The angle 6 is the local inclination of the body 
surface.) If the left-hand and right-hand terms of Eq. (1) are 
divided, respectively, by C and by its equivalent for blunt 


Pmax 


noses C, sin? 6 there results 


mar max’ 


C,/C> 


= sin? 6/sin? 6 nar (2) 


This equation will be referred to as the generalized- Newtonian 
theory. 

For bnaz = 90° the utility of the generalized- Newtonian theory 
is obviously well established since in this case it reverts to Lees! 
blunt-nose modification. For 6a: < 90° (pointed-nose bodies) 
the generalized form is equally useful, and in addition to exhibit- 
ing advantages over the Newtonian theory (C, = 2 sin? 6), it 
tends to unify the results for body shapes in general at hyper- 
sonic speeds. To illustrate, Fig. 1 compares exact solutions 
(including rotational effects) for several ogives with the general- 
ized-Newtonian theory and the Newtonian theory. The exact 
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solutions for values of the hypersonic similarity parameter K = 
M(d/1) from 0.5 to 2 were obtained by Rossow.? For ogives the 
ratio sin? 6/sin? dn.az equals [1 — (x//)]? and is independent of the 
fineness ratio //d of the ogive; here x// is the fraction of the axial 
length of the ogive. The generalized-Newtonian theory for 
ogives is thus 


C,/C> [1 — (x/I)]? (3) 


and, as shown in Fig. 1, its implied independence of K is con- 
firmed by the exact solutions with which it is in excellent agree- 
ment, except near x// = 1 where good agreement requires K >1. 
The hatched band indicates the Newtonian predictions for the 
same range of A covered by the exact solutions; the Newtonian 
predictions are markedly affected by the value of K, and they 
do not approach the accuracy of the generalized-Newtonian 
theory except for K > 1. Note that in predicting C, by use of 
the generalized-Newtonian theory, the exact value of Cp, is 


maz 


to be used. 

Fig. 2 illustrates more graphically the unifying concept em- 
bodied in generalized- Newtonian theory, in that this compilation 
includes both two-dimensional shapes (filled symbols) and axi- 
symmetric shapes (open symbols) for both attached and de- 
tached shocks. The ratio of specific heats y for these data is 7/5 
unless otherwise indicated. The general correlation of the 
data is good and, with few exceptions, there is close agreement 
with the generalized-Newtonian theory. The significant ex- 
ceptions appear to be associated either with low free-stream 
Mach Numbers, where the theory is not expected to apply, or 
with blunt two-dimensional surfaces. With regard to the latter, 
there is some indication that for 17 > 1 the data would parallel 
the generalized-Newtonian theory except for values of sin? 6/sin? 
bmazx near 1. The values of 6,4, used for the cone and flat plate 
(or wedge) data are arbitrary selections and have no special 
significance. 

The generalized-Newtonian theory does not experience the 
significant loss in accuracy with increasing y that occurs with 
Newtonian theory, particularly in application to two-dimensional 
surfaces. In fact, by use of hypersonic small disturbance theory 
for >1, ie. Cy, = (vy + 1367; an analogy may be drawn with 
generalized- Newtonian theory to show that generalized-New- 
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Fic. 2. Correlation for two-dimensional and axisymmetric 
shapes by generalized-Newtonian theory. 


tonian theory is independent of y. Two sets of data for y = 
5/3 are included in Fig. 2 for illustration. 
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The Matching of the Viscid and Inviscid 
Regions for the Stagnation Magnetic Flow 


Paul S. Lykoudis 
Associate Professor, School of Aeronautical Engineering, 
Purdue University, Lafayette, Ind. 


December 12, 1958 


A CONSIDERABLE NUMBER OF PAPERS have been published re- 

cently on the stagnation flow around a blunt nose at hyper- 
sonic flight in the presence of a magnetic field fixed relative to 
the moving body and acting perpendicularly to the streamlines 
of an ionized fluid. 
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In the studies of references 1-3, no attempt is made to calcu- 
late the distortion of the bow shock wave in front of the blunt 
body and, therefore, to relate the inviscid magnetic flow to the 
magnetic boundary layer. This incongruity was pointed out 
in reference 4, whereas recently in reference 5 a solution to the 
magnetic inviscid case was offered. In the light of this last 
work, some comments are pertinent in matching the magnetic 
viscid and inviscid regions. 

In reference 3, the general problem of the existence of simi- 
larity solutions for wedge types of flow was discussed for com- 
pressible ionized fluids in the presence of a magnetic field acting 
perpendicularly to the direction of flow. With the usual bound- 
ary-layer notation, the following are the equations of momentum 
and energy conservation: 


= + 3) ~ (1) 

S’ + Prf S'=0 (2) 

where = (3) 
For similar solutions it is required that 

(4) 


For the stagnation case m = 1 and B may be assumed to be con- 
stant along the direction of flow XY. 
From Eq. (1) one may deduce that 


where C is the inviscid constant and: 
Z = (¢/2m) + V1 + (¢/2m)? (6) 


It has also been found that the solution of Eqs. (1) and (2) may 
be reduced to the solution of the compressible nonmagnetic 
case (subscript zero) by making use of the transformation: 


= fo’ — (1 — (7) 
For cool walls this reduces approximately to: 
= (8) 
The energy equation becomes: 
+ PrZfyS' = 0 (9) 


and hence it is obvious that for the stagnation flow case we have: 
q= Qmag. /Qnonmag. (10) 


For m = 1, Cis the gradient of the velocity at the body and may 
be calculated from reference 6 or, more simply, from the following 
approximation: 
C = V po/pol2 — (px/po)] (11) 
The quantity c is the radius of curvature of the shock. In 
order to test the accuracy of Eq. (10) with the findings of refer- 
ence 5, it is appropriate to calculate the magnetic parameter 
given by Eq. (3) for an average magnetic field over the distance 
between the surface of the body and the position of the shock 
wave. Since in reference 5 the assumption is made that B is in- 
versely proportional to the cube of the radius of curvature, we 
calculate easily that: 

Bavg. = (A + 2)/[2(A + 1)?]Bo (12) 
where A is the nondimensional shock stand-off distance and 
By is the magnetic intensity at the stagnation point. 

The velocity gradient at the wall is given in reference 5, as 
follows: 


Cras. = (uo/c)Vs (13) 


The function V; is tabulated for pop/p2 = 11 and different values 
of the quantity Q, which is related to the present parameter ¢ as 
follows: 


Po Vs poVs 
po OA + 1)? (A + 2)?L + 1)? (A + 2)? 
The ratio of Cmag. computed from Eqs. (5) and (13) for the 


+> 1] (14) 
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| TABLE 1 
| 0 0 0.067 0.524 1.00 1.00 
1.00 0.93 0.088 0.320 1.04 1.02 
2.00 2.12 0.131 0.198 1.05 1.02 
| 3.00 4.05 0.233 0.120 1.02 1.01 
3.40 5.50 0.299 0.088 1.05 1.02 
| 38.80 8.95 0.477 0.057 1.02 1.01 
| 4.00 15.80 0.754 0.034 (0.98 0.99 


values po/po = ll and Rey = 1 is equal to: 


= Venonmag./V Vel Vs + (pe/po) (A + 1)2(4 + (15) 


This ratio is given in Table 1 for different values of the param- 
eters, Q, and the corresponding values of A, Vs, and ¢. 

From Table 1, it is seen that the agreement between the two 
methods of calculation of the velocity gradient at the wall is ex- 
cellent for all values of Q. The following conclusions may now 
be drawn: 

Eq. (5) is a good approximation for the calculation of the ve- 
locity gradient at the wall in stagnation flow, provided that an 
average magnetic field intensity is chosen, as indicated by Eq. 
(12). Inasmuch as \ = 1, Eq. 5 suggests a simple relation yield- 
ing the ratio Vs/Vsnm as a function of (p~/po), A and Q. 

Vs 1 po 

=-—- — (A + 1)?(A + 2)?0+ 


9 
snm po 


2 9 )2, 2 
This is compatible with the findings of reference 3. 

If use is made of the approximation of Eq. (11), then, taking 
into account Eqs. (3), (10), and (11), we can calculate the intensity 
of the magnetic field in Gausses at the stagnation point as a 
function of the flow parameters and the nondimensional heat 
transfer parameter @ as follows: 


B, = 10,000 + Ju +1 Va + 
Now NoDN {(4/2) + 1} 


(17) 
where “#. is in ft./sec. and oD is in mhos. D is the diameter of 
curvature at the nose of the body. 

In the light of Eq. (17), it is of interest to investigate the rela- 
tive effect in the calculation of the magnetic field intensity of the 
following two causes: (a) alteration of the inviscid pressure dis- 
tribution because of the existence of the magnetic field; and 
(b) weakening of the magnetic field intensity in the radial direc- 
tion away from the wall. 

The following three factors may be formed: 


m= V/VA+1, (44 1/{(4/2) + 
Az = (18) 


These factors represent the ratio of the exact value for the mag- 
netic field intensity, divided by the magnetic field calculated 
by: (1) taking into account (b) but not (a); (2) taking into 
account (a) but not (b); and (3) neglecting the influence of both 
(a) and (b). 

These values are tabulated in Table 2. 

It is seen that Case 1 overestimates the magnetic field intensity 


TABLE 2 
e As 
0 0.97 1.10 1.07 
1.00 0.96 1.13 1.08 
2.00 0.94 1.20 1.13 
3.00 0.90 1.36 1.22 
3.40 0.88 1.47 1.29 
3.80 0.82 1.76 1.45 
4.00 0.76 2.24 1.69 
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to produce a given reduction in the heat-transfer rate at the wall, 
whereas the contrary is true for Cases 2 and 3. It is furthermore 
seen that the two causes (a) and (b) above have a reverse effect 
in the calculation of By and that cause (b) is by far more in- 
fluencing than cause (a). This suggests that the weakening of 
the magnetic field (and, therefore, of the resultant ponderomotive 
force) with the inverse third power of the distance from the pole, 
is much more important than the alteration suffered in the pres- 
sures prevailing around the body. 
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Note on the Determination of Higher Modes 
of Vibration by the Stodola or 
Matrix-Iteration Method 


Henry E. Fettis 
Applied Mathematics Research Branch, Aeronautical Research 
Laboratory, WADC, ARDC, Wright-Patterson AFB, Ohio 


December 12, 1958 


A COMMONLY USED METHOD for finding the natural frequencies 

and the associated deflections of beams or other elastic 
structures is the process of successive approximations to the de- 
flection, originally devised by Stodola for computing vibration 
modes of shafts and since adapted to similar problems by numer- 
ous later investigators. 

Basically, all forms of the method are the same. An initial 
assumption for the deflection is made, and a new deflection is 
computed either by numerical integration or from a matrix of 
influence coefficients. Repetition of the process results in con- 
vergence to the fundamental mode—i.e., to the mode having 
the lowest natural frequency. The convergence is dependent 
on the fact that any assumed deflection is a linear combination of 
all of the normal modes, and that, by the iteration process, the 
higher modes are de-emphasized in favor of the dominant one. 
In order to obtain convergence to a nonfundamental mode, it is 
customary to suppress, or ‘“‘sweep-out’’ the lower mode (or 
modes) by means of an orthogonalization process. This, as it is 
well known, requires that the associated deflections, as well as 
the lower frequencies be known to a high degree of accuracy. 
It may also require—in cases where the matrix is nonsymmetric— 
the computation not only of the so-called direct modes but also of 
the associated or conjugate modes. 

The purpose of this note is to show how the lower modes can 
be eliminated if merely the frequencies of these modes are known, 
and, in fact, even when they are not known to a high degree of 
precision. Thus, the numerical complications resulting from the 
enforcement of orthogonality conditions are completely avoided, 
and the process can be applied without modification to either 
symmetric or nonsymmetric matrices. 

Let us assume that we start with an assumed deflection, y, and 
compute from it a new deflection y™. Because of the complete- 
ness property of the normal modes y;, the assumed deflection can 
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be expressed in term of the normal modes: 
y = ay: + aye +... + (1) 
The new deflection is then 
y) = ary + + (2) 
where the \’s are the reciprocals of the squares of the natural 
frequencies taken in ascending order. If A; were known exactly, 
we could form from these two vectors, a new one, Z, such that 
z= — = Ay — Ay (3) 
Clearly, 
z= ary A2)Y2 + + — xX, Y, (4) 
and would therefore be free of the y; component. We could now 
use Z as an iterating vector, and obtain the succession of vectors 
2) = — Ao +... + — An Wn } 
2°?) = — +... + Gndn2Qr — (5) 
= — Az )Y2 +...+ — Wr! 
which would clearly converge to the second mode. In practice, 
an exact value of \; can never be obtained, so that it is necessary 
to repeat the original process of removing the contribution due 
to y; after several iterations with z. A better insight into the 
exact nature of the process is obtained by considering the asso- 
ciated matrices involved. Thus, Eq. (3) is equivalent to operating 
with the vector y on the matrix. 
(6) 
which has the eigenvalues: 
O, — As, At — An 
However, these are arranged in the opposite order from that which 
would normally be desired. On the other hand, the matrix 
B, = A(Al — A) (7) 
has the eigenvalues 
and, in general, the matrix 
B, = A*(AI — A) (8) 
will have eigenvalues 
0, — Ae), — An) 
and it is readily shown that it is always possible to determine a 
value for k, such that 
— > Ast(Ar — As) (9) 
so that second eigenvalue of the original matrix will be brought 
into the dominant position. However, the maximum value of 
k is restricted by the accuracy with which \; is known. This is 
due to the fact that the eigenvalues of B, are, in reality, 
where }; is the known approximation to \;, and it is now evident 
that k must also be so restricted that 
— < — (10) 
in order that the dominant eigenvalue of the original matrix 
remain suppressed. Actually, this is not a severe restriction if 
»; and A» are reasonably well separated. For example, with 
A = 3, A» = 2, = 2.9, Eq. (10) is satisfied for k < 5. Nor- 
mally, a much more accurate value for \; would be known and, 
therefore, a much larger value of k could be admitted. 
Continuation of the process, with \; and \» known, requires 
forming the vectors 
— Ay! 
(11) 
w = dow — Az) 
In theory, at least, w is of this form 


w= dz) Qe As Ys +...-4+ 
alm An) (re An Yn, (12) 


so that the sequence 
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= (13) 
would lead to convergence to d;, if 4; and A, were exact values. 
However, as before, the occasional removal of the contributions 
from y; and y2 is necessary due to inherent numerical inaccuracies 
in the values of \; and do. 

The above process may be continued until all the eigenvalues 
have been found. 

The advantages of the present scheme for either hand or ma- 
chine computation are obvious. Of primary importance is the 
fact that the vectors themselves never enter directly, and also 
that a high degree of precision in the eigenvalues is not even 
necessary. In fact, for matrices of relative low order (say, eight 
or less), it appears that the most efficient procedure would be to 
obtain relatively crude approximations to all of the eigenvalues, 
and then to use these values as a means of systematically im- 
proving the approximate ones. 


w) = Aw 


The Prospects for Magneto-Aerodynamics— 
Correction and Addition* 


E. L. Resler, Jr., and W. R. Sears 

Professor of Aeronautical Engineering, and Director, Graduate 
School of Aeronautical Engineering, Respectively, Cornell 
University, Ithaca, N.Y. 

December 29, 1958 


NE OF THE MAJOR PURPOSES of our paper! was to try to pre- 
dict the probable orders of magnitude of the major electrical 


* This work was carried out at the Graduate School of Aeronautical Engi- 
neering under the sponsorship of the Mechanics Branch, Office of Naval Re- 
search. Reproduction in whole or in part is permitted for any purpose of the 
United States Government. 
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Fic. 1. Electrical conductivities of air and air seeded with 
0.1 per cent potassium behind normal shock wave, as functions 
of the shock Mach Number. Pressure before shock is 1 mm. 
Hg. Data for air are from reference 2, and data for seeded air 
have been computed from these as described herein. 


effects. It is particularly unfortunate, therefore, to find that we 
introduced an extraneous factor of 1/100 in plotting the diagram 
(Fig. 2) in which the electric body force is compared with the 
typical aerodynamic force. All of the values of the ordinate in 
Fig. 2 must be multiplied by 100; thus the electrical force may 
actually become large compared to the dynamic-pressure effect 
if the assumed numerical values of conductivity, field strength, 
etc., can actually be provided. The error was called to our at- 
tention very kindly by James C. Williams, III, of the University of 
Southern California. 

Since this correction makes the result so much more optimistic, 
it might be well to remind the reader that the calculation is in- 
tended to suggest orders of magnitude only. The use of the 
dimension L (e.g., Z = 1 m. in our diagram) implies that the 
body force can be maintained constant throughout a region of 
space having this typical dimension. 

We also take this opportunity to call our readers’ attention to 
important new data on the conductivity of air, ¢. These data, 
presented in Fig. 1 herein, have been obtained by Dr. S. C. Lin? 
in shock-tube experiments. They are more reliable than the data 
presented in Fig. 1 of our original paper, which were based on 
computations using then available information on cross sections 
of air constituents. As indicated here, Lin’s data are for an 
initial pressure of 1 mm., so that they correspond to an altitude 
of flight of about 150,000 ft. At comparable conditions his data 
show considerably higher conductivities than we predicted, and 
we suggest that these more optimistic values be used in place of 
our original Fig. 1 for estimates of air conductivities. 

We have used Lin’s data to estimate a revised curve of o for 
seeded air; the result is also shown in Fig. 1 of this note. To 
make this estimate, we assumed that Lin’s data constitute new 
and better information on constituent cross sections, notably the 
cross section of oxygen atoms, which are most important in de- 
termining the conductivity in this range of temperatures. We 
have therefore recomputed the seeded-air data of reference 1 
using the cross sections implied by the experimental results. The 
conductivities predicted are now several times as great as those 
we presented in our paper. 

The effects of increased values of o are significant in all the 
numerical predictions we made in our paper. The “prospects,” 
in short, seem to be brighter. 
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An Approximate Solution for Transonic Flow 
in Cascadest 


Albert D. Wood* and Joseph H. Clarke 


Instructor, and Associate Professor, Respectively, 
Division of Engineering, Brown University, Providence, R.I. 


December 19, 1958 


gpm STEADY, ISENTROPIC, isoenergetic, transonic flow through 
an unstaggered, nonlifting, two-dimensional cascade of 
slender airfoils is governed by the familiar small-disturbance 
equation 

where ¢ is the perturbation velocity potential, y is the ratio of 
specific heats, M,, is the incident Mach Number, and U’,, is the 


+ This research was supported by the U.S. Air Force under Contract No. 
AF 18(600)-664 monitored by the AF Office of Scientific Research of the 
Air Research and Development Command. 

* Now with the Research and Advanced Development Division, Avco 


Manufacturing Corporation, Wilmington, Mass. 
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Fic. 1. Minimum surface pressure coefficient vs. incident 
Mach Number for an unstaggered, nonlifting cascade of circular- 
are airfoils. Cascade solidity = 2/3, and airfoil thickness ratio 
= 0.1 


incident stream velocity. Since the leading coefficient can be 
written (1 — 1/2), the nonlinear term evidently accounts for local 
deviations of M? from the incident value. Eg. (1) assumes a 
form affinely related to Poisson’s equation 


(1 — + = +1)/U 0) Ma? = f(x, 9) (2) 


when f(x, y) is regarded as given, and may be recast as an in- 
tegral equation through use of either Green’s theorem or double 
exponential Fourier transforms. Introduction of the boundary 
conditions at infinity and on the airfoil surfaces then gives, for 
M.. <1, the result 


ihe f f dn X 


2) (0/dx) In [(x — + (1 — — (8) 


The subscript L denotes the usual linear subsonic solution of the 
homogeneous equation. The second term is a summation of the 
influence of sources distributed over the entire flow field, whose 
unknown strength f(£, 7) accounts for the nonlinear effects which 
are important when local Mach Numbers are near 1; these 
sources do not influence the surface boundary conditions. 

Because small-perturbation transonic flows exhibit the property 
0/dy < 0/dx, the approximation 0f/dy = 0 in the integrand is 
suggested. For the problem considered, this evidently corre- 
sponds to a one-dimensional flow through the cascade. If the 
inner integration in Eq. (3) is carried out and f(£) is replaced by 
its definition, an exact differential appears; for points within a 
passage, the result is simply 


col, 9) = ¥) — (My +1) — (4) 


where the subscript 1 refers to the one-dimensional solution, and 
the usual approximation for the pressure coefficient cp has been 
used. 
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In Fig. 1 the variation with 1/,, of the pressure coefficient at 
the point of maximum thickness, as predicted by the present 
method, is compared with experiment and with other theoretical 
predictions for the case of a cascade with solidity 2/3 composed 
of circular are airfoils with thickness ratio 0.1. Shown in the 
figure are the results given by the subsonic approximate method 
of Spreiter,' linearized subsonic theory, and simple one-dimen- 
sional flow. Also shown is the critical pressure coefficient, which 
indicates the pressure for a local Mach Number of unity. The last 
two relations named have been subjected, for consistency, to the 
transonic small-disturbance approximations, and this one- 
dimensional relation is the one used in Eq. (4). The presented 
adaptation of Spreiter’s method, originally developed for an 
isolated airfoil, involves a further approximation for the influence 
of the additional airfoils, but its effect on the final result is negli- 
gible. It is seen that the present solution and Spreiter’s solution 
coalesce with the linear result at low 1/,, as they should, and 
terminate, along with the one-dimensional solution, with vertical 
slopes at their respective cutoff points. The chordwise variation 
of c, for various M,, and other relevant details are presented in 
reference 2. 


The experimental results were obtained in the Brown Univer- 
sity transonic tunnel by simulating the cascade as the images of a 
single airfoil placed between two solid walls; these walls, ideally 
parallel, were, in fact, slightly diverged to allow for wall bound- 
ary-layer growth.? The airfoil boundary layer was tripped to 
minimize errors due to separation. The data were reproduced 
twice within a deviation too small to distinguish the three sets of 
points on the figure. The principal error is due to nonsimulation 
of the theoretical model because of boundary-layer displacement 
effects on airfoil and walls. It is difficult to estimate these effects 
near the critical and choked conditions where their influence on 
the results is probably the greatest, but the corrections would 
shift each point somewhat down and to the left. The usual 
transonic tunnel interference errors are absent in cascade flow 
because the tunnel walls are proper boundaries. Thus, the data 
presented are judged adequate, with qualifications in the region 
where the flow is virtually choked. 


The solution (4) terminates in the figure at the incident Mach 
Number for which the cascade is choked according to the second 
term. This Mach Number exceeds, as it should, the critical 
Mach Number predicted. The mixed flow given by the method 
corresponds physically to symmetric, shock-free recompression. 
The experimental chordwise pressure distributions? indicate that 
the asymmetry introduced by the actual recompression-with- 
shock is small for the geometry considered. Spreiter’s solution, 
on the other hand, terminates at the critical incident Mach 
Number. The present method and Spreiter’s method both 
follow the experimental trend rather well and it is, therefore, of 
interest to remark that the latter method is also interpretable in 
terms of a related one-dimensional argument, which, however, 
is applied instead to the coefficient of ¢,, in Eq. (1). The fact 
that the present prediction is closer to the experimental curve is 
not inconsistent with Spreiter’s finding! that his solution tends 
to overestimate the value of (—Cy min). 


It appears that the simple method presented yields an accu- 
rate solution for the cascade flow considered and should be appli- 
cable as well to the throat region of certain effusers and to other 
internal flows. Perhaps it also gives some insight concerning 
the significance of Spreiter’s method, which is based on an ap- 
proximation formally quite different. 
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Determination of Local Skin Friction 
by Means of Surface Total-Pressure Probes 


James A. Martin 

Aeronautical Research Engineer, High-Speed Flight Station, 
NASA, Edwards, Calif. 

January 6, 1959 


obs ORIGINAL INVESTIGATION by Preston! and extensions by 
Hsu,? and Fenter and Stalmach* have demonstrated that so- 
called ‘‘surface impact-pressure probes”’ (correctly called ‘‘surface 
total-pressure probes”’ or, simply, ‘‘Preston tubes’’) may be used 
for the determination of local skin friction in turbulent boundary 
layers for isothermal, isobaric conditions. The work of Fenter 
and Stalmach was recently summarized here. The present 
writer believes that related work of J. F. Naleid® should also be 
brought to the attention of the readers of this magazine. 
Naleid has further extended the method to the case of com- 
pressible flow with adverse pressure gradient and zero heat trans- 
fer. The tests performed by Naleid were at a nominal free- 
stream Mach Number of 2 with pressure gradients ranging from 
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a favorable gradient of dp/dx = —0.15 psi/in. to dp/dx = +0.69 
psi/in. The experiment by Naleid utilized the wind tunnel used 
for the experiments reported in reference 3, but with a modi- 
fication incorporated to induce an adverse pressure gradient along 
the test section length. 
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the corresponding number’ being mentioned in the text. The 
= arrangement should be as follows. For books: ' Durand, W. F., 
Aerodynamic Theory, 1st Ed., Vol. 1, p. 23; Julius Springer, 
Berlin, 1934. For magazines: "England, C. R., Crawford, 
=A B., and Mumford, W. W., Some Results of a Study of Ultra- 
Short-Wave Transmission Phenomenon, Proc. IRE, Vol. 20, 
No. 12, pp. 481-482, March, 1933. Please give author, title, 
¢dition, volume, number, page, publisher, and place and date of 
publication as indicated. Omission of one required fact causes 


Si Much extra editorial work and possible inaccuracies. 


ILLUSTRATIONS: Illustrations should accompany manuscripts, 
and each should always be referred to in the text by number. 
Drawings or graphs should not be larger than 12 X 16 inches, 
and must be made with jet black India ink on white paper or 
tracing cloth, the latter being preferred. Do not use typewriter 
for lettering. The smallest lettering on 8 X 10 inch figures should 
be no less than '/, inch high. Cross-section paper (white with 
black lines) may be used, but it should not have more than 4 lines 
per inch. If finer ruled paper is used, the major division lines 
should be drawn in with black ink, omitting the finer divisions, 
In the case of finely ruled paper, only blue-lined paper can be 
accepted. Tracing paper and blueprints are not acceptable, 
Lettering and all markings must be large enough to be readable 
after reduction to single-column width (3*/\,in.). Mail rolled or 
flat; never fold. Drawings that cannot be reproduced (including 
pencil drawings) will be returned to the author for redrawing, 
thus delaying publication of the paper. Photographs should be 
distinct and show clear black and white contrasts. They must 
be on glossy white paper. Avoid round and oval photographs, 


CAPTIONS AND LEGENDS: Legends or captions must accompany 
each drawing or photograph submitted. If written on the draw- 
ing or photograph, they should be placed below and well outside 
the part to be reproduced. Each table should have a caption 
such as Table 1, Table 2, Table 3, etc. Captions should be com- 
plete in themselves so as to make the data intelligible to the 
reader without reference to the text. A duplicate list of captions 
for figures should be included as the last page of the manuscript. 
Use “Fig 1’’ (not Figure 1), ‘Figs. 3 and 4,” etc., in both the text 
and the numbering of illustrations. In the text, “Eq. (1)” or 
“Eqs. (1) and (2)” is used, not ‘‘Equation (1).” In captions, 
legends, and in table headings, write all words in full; do not 
abbreviate, except for ‘‘Fig.” and “Eq.” 


MATHEMATICAL WorRK: Formulas may be typewritten or 
carefully written in pen and ink, the writing to be large enough so 
that it can be marked for the printer. Considerable space for mark- 
ing should be allowed above and below all equations. All compli- 
cated equations should be repeated on separate sheets with plenty 
of space left for marking. The solidus should be used for sim- 
ple fractions appearing within the text. Make all expressions 
clear to the typesetter. Greek letters used in formulas should 
be clearly designated by name in the margin of the manu- 
script. The difference between capital and lower-case letters 
should be clearly distinguished and care taken to avoid confu- 
sion between zero (0) and the letter (0), between the numeral 
(one) and the letter (ell) and the prime (’), between alpha and 
a, kappa and k, u and mu, v and nu, n and eta. All sub- 
scripts and exponents should be clearly distinguished. Avoid 
complicated exponents and subscripts. Dots and bars over 
letters or mathematical expressions should also be avoided. 
When it is necessary to repeat a complicated expression, it should 
be represented by some convenient symbol. 


SyMBOLS AND ABBREVIATIONS: The symbols recommended in 
the American Standard Association ‘“‘Letter Symbols for Aero- 
nautical Sciences,"” ASA Y10.7—1954, should be used wherever 
practicable. All symbols should be clearly written and carefully 
checked. Standard abbreviations should be used, and it should 
be noted that most abbreviations are lower case, such as m.p.h., 
b.m.e.p., i.h.p., b.hp., hp., etc. 
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